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ABSTRACT 


This thesis considers the problem of time-optimal spacecraft slew 
maneuvers. Since the work of Bilimoria and Wie it has been known that the time- 
optimal reorientation of a symmetric rigid body was not the eigenaxis maneuver 
once thought to be correct. Here, this concept is extended to axisymmetric and 
asymmetric rigid body reorientations with idealized independent torque 
generating devices. The premise that the time-optimal maneuver is not, in 
general, an eigenaxis maneuver, is shown to hold for all spacecraft 
configurations. The methodology is then extended to include spacecraft control 
systems employing magnetic torque rods, a combination of pitch bias wheel with 
magnetic torque rods, and finally to control systems employing single gimbal 
control moment gyros. The resulting control solutions, designed within the 
limitations of the actuators, eliminate the requirement to avoid actuator 
singularities. Finally, by employing sampled-state feedback the viability of real¬ 
time optimal closed loop control is demonstrated. 
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I. INTRODUCTION 


A. MOTIVATION 

Modern satellites require more precision and agility than ever before. In 
both military and commercial applications the ability to rapidly maneuver 
represents an increase in mission effectiveness and productivity. Rapid 
retargeting maneuvers translate directly into more time on the intended object 
and more observations per orbit. With commercial Earth-imaging satellites like 
Ikonos and research into agile microsatellites 1 on the rise, the need for time- 
optimal satellite control is greater than ever. Planned future satellite missions in 
support of missile-defense and Earth observation will rely on satellite agility and 
minimum-time maneuvers for mission success. 

B. THE PROBLEM 

Spacecraft time-optimal attitude maneuvers have held the interest of 
engineers and mathematicians for decades. In their paper, “Survey of Time- 
Optimal Attitude Maneuvers,” Scriverner and Thompson provide a summary of 
work in this and related areas 2 . They also present the principal difficulties 
associated with solving the time-optimal attitude maneuver problem. 

Unlike the pointing problem which lends itself well to linearization 
techniques, slew problems, especially large angle slew maneuvers are highly 
nonlinear. Euler’s equations, which represent the rotational motion of a rigid 
body, consist of three, coupled, nonlinear differential equations. 

l/0 X +(l 2 -ly)(0/0y=U, 

/y®y + ( 4 - 4 )®/°z = U 2 
l z a z + (ly-l x )o3 y a x = u 3 

Additionally, the high angular velocities associated with time-optimal 
maneuvering cause significant gyroscopic stiffness through the non-linear terms. 


1 



While these equations completely describe the rotational motion of the 
body with respect to an inertial frame they do not determine the spacecraft’s 
attitude. To describe the attitude of the spacecraft Euler parameters are 
generally used since they provide singularity free kinematics: 
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These non-linear ordinary differential equations have no closed-form solutions 
except for a small number of cases involving simple rotations or simple geometry 
(i.e., torque free motion of an axisymmetric body 3 4 ). 

The time-optimal nature of the problem adds additional difficulties. Major 
problems are encountered in solving the boundary value problem that arises from 
the application of the Minimum Principle. The problem, formulated in Chapter II, 
has no known numerical or analytical solution except when certain restrictions 
are applied. These assumptions have taken the form of specific configurations or 
restricted motions 5 . 


C. HISTORICAL BACKGROUND 

A short summary of the historical work that influenced this research is 
presented here. Beginning with their landmark work in 1993, Bilimoria and Wie 
demonstrated with extensive analytical modeling and numerical analysis that the 
time-optimal maneuver was not the long-assumed eigenaxis maneuver 5 . Their 
methods and results are examined in detail in Chapter III. Bilimoria and Wie later 
extended this work to an axisymmetric body and observed the same 
characteristic precession they had noted earlier. Beyers and Vadali 8 reproduced 
the results of Bilimoria and Wie but focused on developing a control algorithm. 
They used linearization and the switch time optimization (STO) algorithm 
developed by Meier and Bryson 9 to produce an algorithm that could be 
implemented in real-time. 
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In the context of magnetic control, Junkins et al. 10 modeled a single 
controller aligned with the spin axis of a spin-stabilized symmetric spacecraft. 
The time-optimal solutions were found through an interactive graphical 
technique. These results were eventually implemented for open loop control on 
the NOVA-1 spacecraft 11 . 


D. CONCLUSION 

Solving the problem of spacecraft time-optimal control has occupied the 
interest of engineers and mathematicians for years. The problem is simple to 
formulate and yet solutions have been difficult to obtain. In this thesis we extend 
the present body of work to all spacecraft moment of inertia configurations. The 
ideal actuators often studied to simplify the problem formulation are replaced with 
magnetic torque rods and a combination of toque rods and a pitch-bias wheel. 
Finally, the control moment gyro configuration most studied for its singularity 
problems is examined and shown to be singularity free in the time-optimal result. 
This unpredicted benefit of the formulation has the potential to make future work 
in actuator singularity avoidance moot. Finally, computational speeds are shown 
to be such so as to allow time-optimal solutions to the nonlinear plant to be 
generated in real-time. 
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II. OPTIMAL CONTROL PROBLEM FORMULATION 


A. INTRODUCTION 

In this section the optimal control problem formulation is developed. This 
formulation will be referred to throughout the remainder of this work. There are 
numerous excellent books and papers on optimal control theory which explain 
these concepts in much greater detail. Interested readers are referred to the 
references [1,2, 3, & 4] for further information. 


B. METHODOLOGY 

Over the years many different methods of solving optimal control problems 
have been developed. These are broadly grouped into two categories: indirect 
and direct methods. In indirect methods, the necessary conditions for optimality 
are derived from Pontryagin’s Minimum Principle and solved to obtain the optimal 
trajectory. These methods are notoriously labor intensive. In direct methods, the 
optimal control problem is discretized into a parameter optimization problem. 
The resulting nonlinear programming problem can then be solved by standard 
nonlinear programming means. 

In this work we will employ a Legendre Pseudospectral method encoded 
in the reusable software package DIDO 5 . Pseudospectral methods are well 
known in the field of fluid dynamics where they are used to numerically solve 
partial differential equations. Unlike other methods which employ piecewise- 
continuous polynomials Pseudospectral methods are unique in their application 
of global orthogonal polynomials as trial functions. 

In the Legendre Pseudospectral method the time domain of the problem is 
discretized at a special set of Legendre-Gauss-Lobatto (LGL) points. Polynomial 
approximations of the state and control variables are considered where Lagrange 
polynomials are the trial functions and the unknown coefficients are the values of 
the state and control variables at the LGL points (nodes). Using the properties of 
the Lagrange polynomials the nonlinear differential state equations are 
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transformed to nonlinear algebraic equations. These equations are then posed 
as a nonlinear programming problem and a sparse numerical optimizer is used to 
solve the problem. In addition, a relationship between the costate variables and 
the Lagrange multipliers called the Covector Mapping Theorem, allows for 
determination of the costates at the LGL points 6 . This provides numerical 
information that is used to validate the solution’s optimality. Interested readers 
should refer to references 7, 8, 9, and 10 for details regarding Pseudospectral 
Methods. 

Numerical results throughout this work are specified in terms of the 
number of LGL points used to obtain the solution. Initial solutions were generally 
based on 30 LGL points with initial guesses restricted to two vectors representing 
the initial and final conditions. Initial control solution guesses were arbitrary. 
Where greater accuracy was desired the 30 LGL point state and control solution 
was used as a guess for a second solution based on 100 LGL points. 

C. OPTIMAL CONTROL PROBLEM FORMULATION 

In an attempt to develop the notation and methodology that will be used 
throughout this work we consider the following optimal control problem. 
Determine the control function u* (t) and the corresponding state trajectory 

x* (if) that minimize the Bolza cost functional* 


JU(),u()J l ) = E{x(t l ),t l ) + ]F(xl,t),u(t),t)dt 

1 0 

where xe R n and ueR m are subject to the differential constraint 

k=f{x(t),4t)) te[t 0 ,tf] 

and the boundary conditions 

§ o {x{t 0 ),t 0 )=0 e 0 eR p and p<n + 1 
§f (x{t f ),t f ) = 0 e f eR q and q<n +1 

(•): Functional dependence not specified. 


6 



and control inequality constraints of the form 

h(u,t) < 0 heR r 

where dh/du has full rank. 

The Lagrange multiplier theory allows us to adjoin the state equations and 
constraints to the cost functional to form an augmented cost functional as follows 

J = E ( x{ t f ), t f ) + Ve 0 (( x(t 0 ), t 0 )+ v f T e f (( x(t f ), t f ) + 

| { F(x, u,t) + X{t) T {f(^,u,t)-x) + [i{t) T h{u,t)}dt 

h 

Defining the Hamiltonian as, 

H(k,x,u,t) = F(x,u,t) + X T f(x,u,t) 

Pontryagin’s Minimum Principle provides the following necessary conditions for 
u* to be an optimal control. 

FI(V,x*,u*,t) <H(X*,x*,u, t ) ->Hamiltionian Minimization 

X* = -> Adjoint equations 

dx 

do 

X(t 0 ) = -v T — L ^ Initial transversality 
dx 0 

X{t f ) = ^-+v> T ^ L -^> Terminal transversality 
dx f ax f 

Fl\t f ] + — +v T — = 0 -> Hamiltonian Value 
L dt dt 

For the case in which the cost function and constraints are linear in the 
state and control variables, no minimum exists for the Hamiltonian minimization 
unless inequality constraints are imposed on the state and/or control variables. 11 
In the case, where the inequality constraints are linear and placed only on the 
control variables the problem is a special case of linear programming problem 
covered in detail by Bryson^ and others 13 . The important result is the application 
of the Karush-Kuhn-Tucker (KKT) Theorem and complementarity conditions. 

Then, if a minimum exists, it will require the control variable to be at a boundary 
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of the feasible control region. Since the Hamiltonian is subject to an inequality 
constraint on the control variable we apply the Karush-Kuhn-Tucker Theorem 
and form the Lagrangian of the Hamiltonian: 

H = H+[L T h (2.1) 


where p, is a KKT multiplier. Then the KKT Theorem gives the following 
necessary conditions for the optimal trajectory 


dH__dH_ 
du : du 


+ 


dh 

du 


(1 = 0 


( 2 . 2 ) 


In addition, the multiplier-constraint pair must satisfy the complementarity 
conditions of the KKT Theorem which states: 


<0 h,(u,t) = ^(t) 

= 0 if h^(t)<h i (u,t)<h] J (t) 
>0 h i (u,t) = h^ J (t) 

unrestricted h[(t) = h“(t) 


(2.3) 


In subsequent chapters, once the state dynamics of the problem have 
been established, be optimal control problem will be presented in this format. 
The candidate optimal control solution will then be propagated through the state 
dynamics to verify the feasibility of the solution. Verifying that this candidate 
solution is optimal is more difficult. Recall that the Minimum Principle supplies 
necessary conditions not sufficient conditions for optimality. However, since the 
Legendre Pseudospectral Method, through the covector mapping theorem, 
provides costate information at the nodes we can evaluate the necessary 
conditions for optimality. Where available, results will be compared to published 
works. Finally, engineering and physical insight will be used to establish that the 
solutions obtained are optimal. 
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III. IDEALIZIED TORQUE ACTUATOR CONTROL PROBLEM 


A. INTRODUCTION 

In this section, we develop the time-optimal rotational maneuver for a 
rigid-body with ideal, bounded actuators. The idealized actuator will be defined 
by: 

Torque out =[l]u 

where \u\ < 1, is the control vector and / is the identity matrix. In the case of the 

ideal actuator, the control vector is equal to the output torque vector. This 
distinction allows for later definition of the control vector u as, for example, 
magnetic dipole moment (See Chapter IV). In this case the actuator is not 
considered ideal. The closest physical approximation to the idealized actuator is 
a thruster. 

Bilimoria and Wie showed that the time-optimal solution for the 
reorientation of an inertial symmetric body was not necessarily an eigenaxis 
maneuver.i Shen and Tsiotras examined time-optimal reorientations of 
axisymmetric spacecraft using two controls. 2 Finally, Proulx and Ross examined 
the control structure and evaluated time-optimal reorientations of asymmetric 
spacecraft 3 . We begin with a reexamination of the inertial symmetric 
reorientation problem in order to establish the methodology. Then the principles 
are extended to axisymmetric and finally to asymmetric spacecraft in the orbital 
frame. 

B. INERTIALLY SYMMETRIC RIGID-BODIES 

Inertial symmetric rigid-body reorientations represent the simplest problem 
that we can pose. Yet, prior to the work of Bilimoria and Wie the solution was 
misunderstood. The eigenaxis maneuver about the control axis was thought to 
be the time-optimal maneuver for symmetric bodies and near time-optimal for 
other configurations. 4 Through analytical analysis and numerical simulation 
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Bilimoria and Wie demonstrated that this was not true for symmetric bodies. 
They showed that all three control components cannot be simultaneously zero on 
the time-optimal trajectory. 


1. Notation and Transformations 

In order to progress with the problem formulation we establish the 
standard notation and rotation sequences used in the development. The 
spacecraft will have an assumed standard orbit frame defined as: 

Z 0 — Nadir pointing 

X 0 - Velocity vector 

Y 0 - completes the right hand set 

We choose the rotation sequence for the body to orbit frame 
transformation, represented as Euler angles as \\r -^>0 -xp with axes order of 
rotation 3 -^2 ->1. Therefore, \\t is the first angular rotation about the zbody 

axis. The second rotation is about the once displaced ybody axis by an angle 
0 . The final angular rotation cp is about the x-body axis. 

Then, following the notation convention of Kane 5 the following angular 
velocities are defined: 


N-b 

COg 



Angular rate of body with respect to inertial in body frame. (3.1) 


N—O A 


CO 


o 


0 

-®o 

0 


-»Angular rate of orbit with respect to inertial in orbit frame. (3.2) 


o — b 

CDo 


co 1 

C0 2 

® 3 


-»Angular rate of body with respect to orbit in body frame. (3.3) 
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and: 


N 65 b B = %y B ° + %J B (3.4) 

The rotational transformation from the orbit frame to the body frame is 
referred to as the direction cosine matrix (DCM). It is defined from the rotation 
sequence above and represented in terms of the quaternion vector as, 6 

4 - 4 - 4 + 4 2 :{q,q z - q 3 q 4 ) 2{q,q 3 - q 2 q 4 ) 

B C°= 2 (q,q 2 -q 3 q 4 ) q\ - q 2 - q 2 3 + q\ 2 {q 2 q 3 -qq 4 ) (3.5) 

2 (q,q 3 -q 2 q 4 ) 2 (q 2 q 3 -q,q 4 ) q\ -q\ -q \+ q\ 

Therefore, the angular velocity of the body with respect to the Newtonian frame 
may be written as: 

X J = °«B + B C° %° 0 (3.6) 

Then substituting from equation (3.2) we have, 

w o) e b = °65 b B -oJ 0 C ( . 2 (3.7) 

where co 0 is the magnitude of the orbital angular velocity* and C /2 is the second 

column of the DCM. Then by substituting the values established in equation 
(3.5) we obtain the following relationships, 

[<0,1 T 2 {q^q 2 +q 3 q 4 ) " 

w y = ca 2 -o) 0 q 2 2 -q]-q 2 3 +q 2 4 (3.8) 

_o) z J [o) 3 J [ 2(q 2 q 3 +q,q 4 ) _ 


Orbital angular velocity co o is sometimes referred to as “mean motion” for circular or 
elliptical orbits. 
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and, 


co 1 


V 


" 2(g 1 <7 2 +q 3 q 4 ) “ 

co 2 

= 

°>y 

+ ®„ 

2 2 2 2 

q 2 - - Qs + Qa 

®3_ 




_ 2(q 2 q 3 + QiQ 4 ) . 


2. Problem Formulation 

The inertial symmetric body is shown in Figure 1. The state of the system 
can be completely defined by its attitude and the time-rate-of-change of attitude. 



Figure 1 Inertial Symmetric Body 

For mathematical simplicity in extension to future work, the attitude is 
represented as a quaternion. Quaternions have none of the inherent singularities 
that are well known in other representations of spacecraft attitude. Additionally, 
they are computationally more suited for on-board real-time processing since 
there are no trigonometric functions to be evaluated in the quaternion kinematics 
equations. 

Euler’s rotation theorem states that a rigid body can be changed from any 
given initial orientation to an arbitrary final orientation by a single rotation about 
an axis that is fixed to the body and stationary in an inertial reference frame. 7 
This axis which remains unchanged in both the body and reference frames is 
called the eigenaxis or Euler Axis. 8 
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Defining the eigenaxis in the body frame as: 

e = ej + ej + e 3 k 

allows us to define the quaternion vector as follows: 

9i = eiSin(%) 

q 2 = e 2 sin|^2 J 

^3 = e 3 sin (%) 

= cos [%) 

where <|) is the rotation angle about the eigenaxis. The state of the spacecraft is 
then represented by:t 



to 


The state dynamics include both quaternion kinematics and rotational 
dynamics. The quaternion kinematics are well known and are repeated here for 
completeness. 


(3.10) 


(3.11) 


q = with, ^ = 


0 

CO 

3 

—CO 2 

co 1 

-co 3 

0 

0^ 

CO, 

c 

co 2 

— 

0 

CO, 

—co 1 

—C0 2 

-co 3 

0 


(3.12) 


Termwise, we have, 


t We have adopted the convention that q 4 is the scalar quantity of the quaternion vector. 
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4 =^[®iQ 4 -«bQ3+c%Q 2 ] 
4 =^KQ 3 + C ^ C ?4-“3 Qi] 

4 =^[-®iQ 2 +“2^1 +“ 3 q 4 ] 

4=^Hi4-^4-^4] 


(3.13) 


In order to maintain consistent notation throughout this work we have used the 
subscripts indicating that the angular velocities are of the body with respect to the 
orbit in the body frame. For this example, we consider the spacecraft located in 
inertial space. Mathematically, orbital angular velocity is zero and we have, 


CO. 


CO 

1 


X 

co 2 

= 

°>y 

® 3 _ 


“z_ 


This distinction will become important for the later case where the spacecraft is in 
orbit about the Earth centered inertial frame. 

The rotational dynamics of a rigid body are obtained by equating the 
applied torque about the center of mass to the time rate of change of the angular 
momentum. 9 This well-known fact is expressed as: 


where H is the angular momentum vector of the rigid-body about its center of 
mass with respect to an inertial frame and M is the external moment acting on 
the body about its center of mass. 

Since this is the logical extension of Newton’s 2 nd law it follows that the 
time-derivative of angular momentum must be with respect to an inertial frame. 
The angular momentum is the product of the moment of inertia and the angular 
velocity, 

H = ia 3 
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This allows us to express Euler’s rotational equations of motion for a rigid-body in 
matrix notation. 

/co + cox /co = M ext 

If we express the moment of inertia and angular velocity in the principal axis 
frame as: 

"/* o o“ 

1 = 0 /, o 

o o / z 

co =co 1 j 6 1 + co>4 + g> 3 4 

Then Euler’s equations can be expanded to: 

^=l x (o x +{l z -l y ) co y co z 

/W 2 =/ y c6 y + (/ x -/ z )co x co z (3.14) 

M 3 =//6 z + (/ y -/ x )o) x co y 

These equations are well known and generally referred to as Euler’s Moment 
Equations. 10 For the case of a symmetric body the gyroscopic terms in equations 
(3.14) are zero and the rotational dynamics can be written as: 

. M . M . M, 

C0x=7 1 c °y = -r 

x y 1 z 

This reduces Euler’s Moment Equations from three-coupled non-linear ordinary 
differential equations to three uncoupled linear ordinary differential equations. 
Then normalizing the moment by the inertia, without any loss of generality we 
have, 

c b x = c6 r = M 2 o6_ = M 3 (3.15) 

Taken together, equations (3.13) and (3.15) form the system dynamic 
constraints. The normalized external torque,M,., is chosen as the control 
parameter and it is limited somewhat arbitrarily to unity. 
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2. Time Optimal Maneuvers 

With the dynamic constraints defined and presuming initial conditions and 
desired final conditions are known, we can formally state the time-optimal control 
problem for the rest-to-rest maneuver as follows: 

Minimize J(x{),u{),t f ) =t f -t 0 

s.t. q = — £lq (3.16) 

co = u, |l/| < 1 

In order to allow a comparison with published results, we choose our initial and 
final conditions as: 


X- 

--[Qi 

9 2 9 3 

9 4 

co 1 

co 2 C0 3 ] 7 


*o 

= [o 

0 0 1 

0 

0 

or 

(3.17) 

Xf 

= [0 

0 1 0 

0 

0 

or 



Referring to equations (3.10) and (3.11) these conditions represent a rest-to-rest 
rotation maneuver of 180 degrees about the z-body axis. 

As our first step in determining the control u* that will drive the dynamic 

system we write the Hamiltonian in standard form. Since we have written the 
cost functional in the Mayer form it will not appear in the Hamiltonian. Thus the 
Hamiltonian will be a linear combination of the state dynamics and take the form: 

H(X,x,u,t) = X T f{x ,u) (3.18) 

Substituting equations (3.13) and (3.15) into the Hamiltonian equation (3.18) we 
have: 

X X 

H = -y(“ i 9 4 - “ 2 9 3 +“392) + y(“i9 b + “>9 4 -® a 9i) + 

X X 

y(-® 1 9 2 +“29i+“ 3 94)+y(-“i9i-“2^-“ 3 93)+ (3-19) 

\o x 9 + X ay u 2 + A , m u 3 
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where the subscripts on the Lagrange multipliers have been selected to aid in 
bookkeeping. 

In order to minimize the Hamiltonian we form the Lagrangian of the 
Hamiltonian by adjoining the constraint equations to the Hamiltonian in the form: 

H = H+[L T h (3.20) 

where |i,are the KKT multipliers and h(u{t),t) is the control constraint function in 

the standard form. On substituting the Hamiltonian equation (3.19) and the 
control constraint equation as defined above into equation (3.20) the inertial 
symmetric problem has necessary conditions, 

h+V = 0 (3.21) 


with the complemetarity condition, 


<0 


Ui = -1 

= 0 

if 

-1 <Uj< 1 

>0 


u. = 1 


(3.22) 


The quantity dH/du is called the “switching function” in the literature. The 

case when the switching function equals zero for a non-zero period of time was 
rigorously examined by Bilimoria and Wie and shown not to be time optimal. 
Thus we are left with a switching function that determines when the optimal 
control u* will switch between its extreme values. For this reason the control 
profile is called bang-bang.n 

In order to validate their results, Bilimoria and Wie used a multiple¬ 
shooting algorithm to solve the two-point boundary value problem resulting from 
the state and adjoint equations. This, in conjunction with the state and adjoint 
equations allowed them to determine and evaluate candidate optimal control 
solutions and the resulting trajectories. 12 Here we employ a pseudospectral 
approximation to arrive at the optimal solution. 
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A careful comparison of previous results with the ones obtained here 
(Figure 2 and Figure 3) reveals that for this case there are at least two and 
potentially four equally optimal solutions. For the symmetric spacecraft the 
precession that is characteristic of the optimal solution may proceed in either 
direction for a 180 degree reorientation. A simple transformation of the controls 
as follows: 

Current Previous 
u 1 -> (—1) x u 2 
u 2 —> q 
u 3 -^u 3 

will transform the current quaternion and angular rate histories figure 2 and 
Figure 3) to an exact match with those of Bilimoria and Wie. These transformed 
results are shown in Figure 4 and Figure 5. For the 180 degree maneuver under 
consideration one could argue that the entire maneuver could proceed in the 
opposite direction and still result in the same cost and hence be equally optimal. 
This, while true, is precluded by the quaternion definition convention employed. 

Therefore, the results obtained correspond to the published results and 
clearly show that the solution is not an eigenaxis maneuver. This is evident from 
both the non-zero quaternion histories of q, and q 2 (Figure 2) and the non-zero 
angular velocities about the x and y axes (Figure 3). 
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Angular Rates History 



Figure 3 Inertial Symmetric Angular Rate Solution 
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Figure 5 Inertial Symmetric Angular Rate Transformed Solution 
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The optimal control solution is shown in Figure 6. This solution matches 
the previous work of Bilimoria and Wie and exhibits the five control switch, bang- 
bang structure that they observed. 



Figure 6 Inertial Symmetric Optimal Control Solution 


The control vector obtained is verified as a feasible control through propagation 
of the state dynamics. The initial conditions and control solution are used as 
input to a MATLAB® ODE45 propagation subroutine which uses an explicit one- 
step Runge-Kutta medium order (4 th - to 5 th -order) solver^ to verify that the 
control solution drives the system from the given initial condition to the desired 
final condition. A linear interpolation was used to approximate the control values 
between LGL points. Propagation results are shown in Figure 7 and Figure 8. 
The original solution obtained is shown in solid lines overlaid with the propagated 
states shown as V marks below. 
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It is easy to see that not only does the dynamic system propagate to the desired 
end state but that the pseudospectral approximation of the states closely 
matches the propagated results. 



Figure 8 Angular Rate Propagated Solution 
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Next we examine the necessary conditions for optimality. Recall that 
equation (3.21) and the complementarity conditions of equation (3.22) define the 
switching structure of the control vector and define a relationship between the 
costate dynamics and KKT multipliers. An inspection of the switching functions 
and their relationship to the control behavior verifies that the control-constraint 
pair meet the KKT conditions. Switching functions for each axis are shown 
(Figure 9, Figure 10, and Figure 11). 



Figure 9 Inertial Symmetric Spacecraft Control and Switching Function 

About X-body Axis 
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Figure 10 Inertial Symmetric Spacecraft Control and Switching Function 

About Y-body Axis 



Figure 11 Inertial Symmetric Spacecraft Control and Switching Function 

About Z-body Axis 
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Turning to the Hamiltonian equation (3.19) there is no explicit dependence 
on time and therefore the Hamiltonian will have a constant value along the 
optimal trajectory. 



Recall the Hamiltonian value condition allows us to determine the final value of 
the Hamiltonian from the end-point Lagrangian. The end manifold, in standard 
form, is given by: 


e{x f ,t f ) 


% 

q 3f -i 

% 

CO. 

7 

co„ 


CO 


3, 


(3.23) 


Then the end-point Lagrangian, previously defined as: 

E = E+v T e (3.24) 

allows us to determine the final value of the Hamiltonian as: 


H [ t f] = _1 


(3.25) 


Thus, the value of the Hamiltonian will be -1 at all times along the optimal 
trajectory. The resultant Hamiltonian, shown in Figure 12, meets the necessary 
conditions with some small numerical variation. 
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Figure 12 Inertial Symmetric Problem Hamiltonian Evolution 

The Adjoint equations can be formed by differentiation of the Hamiltonian. 
However, since the state variables are specified at both the initial and final 
conditions the adjoint variables will be free or unspecified at both initial and final 
conditions. Therefore, the adjoint equations and terminal transversality of the 
adjoint variables provide no new information which will aid in our solution to the 
problem. 


3. Numerical Considerations and Notes 

We have shown that the optimal control solution is feasible and meets the 
necessary conditions derived from Pontryagin’s Minimum Principle. Additionally, 
we have shown that the results obtained closely match those previously obtained 
by Bilimoria and Wie. Next we compare the time required for optimal 
reorientation with that required for a comparable eigenaxis rotation. 
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Symmetric Maneuver Comparison 




45 

90 

135 

180 

□ Optimal 

1.752 

2.428 

2.895 

3.254 

■ Eigenaxis 

1.777 

2.513 

3.078 

3.555 


Table 1 Comparison of symmetric spacecraft reorientation time for time- 

optimal and eigen axis maneuvers 


A comparison of the time required for reorientations is shown in Table 1. 
These results indicate that the time-optimal maneuver is always faster than the 
optimal eigenaxis maneuver except as noted by Bilimoria and Wie 14 . The cost 
reduction is shown below in Table 2. It is clear that larger reorientation 
maneuvers represent a larger cost benefit. 

Proper scaling of numerical problems is important to both the accuracy of 
the result and the computation time required in obtaining the result. No 
numerical scaling was employed in this algorithm. The nature of the problem is 
such that for reasonably large rotational maneuvers all numerical values are of 
the same order of magnitude. 
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Percent Reduction 



Table 2 Cost reduction for inertial symmetric time-optimal reorientations 

over eigen axis maneuvers 


4. Conclusions 

The published work by Bilimoria and Wie represents a landmark 
achievement in the application of optimal control theory. Their basic result, that 
the time-optimal maneuver is not, in general, an eigenaxis maneuver will be 
shown to extend to a wide variety of spacecraft moments of inertia and control 
configurations. 


C. AXISYMMETRIC SPACECRAFT REORIENTATIONS 

In this section we will examine the time-optimal reorientation of 
axisymmetric spacecraft. Shen and Tsiotras 15 examined the problem of 
axisymmetric reorientations using two control torques. They used a combination 
of direct and indirect methods to numerically evaluate several representative 
maneuvers 16 . We begin with reorientations about the axis of symmetry before 
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examining the general reorientation case. Finally, we will examine the 
reorientation of the axis of symmetry of a spacecraft spinning about its axis of 
symmetry. For the spinning spacecraft we will examine the case where the spin 
rate is held constant throughout the maneuver. This amounts to the two-control 
reorientation of Shen and Tsiotras. Additionally, we will examine the case where 
the spin rate is allowed to vary during the maneuver from a given initial condition 
to a known final condition. We will show that adding this third control torque 
results in the time-optimal maneuver for a spinning spacecraft. 

1. Problem Formulation 

The inertial axisymmetric body is shown in Figure 13. The axis of 
symmetry is chosen as the z-body axis somewhat arbitrarily though this is not an 


Z-Body 



X-Body Y-Body 

Figure 13 Inertial Axisymmetric Body 

uncommon configuration. As before the state of the spacecraft will be defined 
as: 



31 



The quaternion kinematics remain as previously defined in equation (3.12). 
However, the rotational dynamics in the principal axis frame given in equation 
(3.14) now take the following form: 

/Wi =/ x oi x +(/ z -/ y )co y co z 

M 2 = /y°i y + (/ x - 4 )ff> x co z (3.26) 

M 3 = / z co z 


Note that the gyroscopic term about the zbody axis has been eliminated by the 
equal moments of inertia about the x and ybody axes. Thus, the time rate of 
change of angular rate is given by: 


co =— 1 


M, 

7 


co, 


CO. 


M 2 

'y 

Ms 

/, 


f i -I A 

z y 

v x ) 

r i-I ' 


|co y co z 


to x co z 


V 7 J 


(3.27) 


These, taken together with the quaternion kinematics equations defined earlier 
form the new dynamic constraints on the axisymmetric system. Once again the 
control parameter is chosen as external torque and is limited to unity. 
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2. Time Optimal Maneuvers 

The formal statement of the optimal control problem is as follows: 


Minimize J(x{),u{),t f ) = t f -t 0 
s.t. q=\&q 


“x = —“ 

'/ -/ 

z y 

* / 

1 

X 

V x 


f, . \ 


h-h 

co 

1 

y 

\ y J 

U 3 


CO = — 


Z / 


z 


|u| < 1 



(3.28) 


a. Non-spinning Axisymmetric Reorientations About the 
Axis of Symmetry 

We will begin by considering a family of prolate spacecraft with 
characteristics as given (Table 3). Assuming the spacecraft has uniform mass 
distribution, the moments of inertia about the principal axes are given by: 



l x = / = — mr 2 +—mf 2 
x y 4 12 


(3.29) 


where, t is the length of the axis of symmetry, mis the mass of the spacecraft 
and r is the radius. The effect of increasing length on moment of inertia is shown 
in Table 3. 
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Case 

Mass 

Length 

Radius 

lx 

iy 

Iz 

1 

1 

1 

1 

0.33333 

0.33333 

0.5 

II 

1 

5 

1 

2.33333 

2.33333 

0.5 

III 

1 

10 

1 

8.58333 

8.58333 

0.5 

IV 

1 

20 

1 

33.5833 

33.5833 

0.5 

V 

1 

50 

1 

208.583 

208.583 

0.5 

VII 

1 

100 

1 

833.583 

833.583 

0.5 


Table 3 Prolate Spacecraft Characteristics 


The initial and final conditions will define a rest-to-rest maneuver in inertial space 
representative of a 135 degree rotation about the axis of symmetry. The initial 
condition will be defined as nadir pointing therefore we have the following: 



q 2 

q 3 q 4 

(0, 

C0 2 

0) 3 ] 7 

*0 = 

[0 

0 

0 1 0 

0 

of 



0 

0 

. f 135° 

sin - 

2 

X 

f 

COS 

135° ^ 
2 


T 


0 0 0 


Referring to equation (3.18) we can write the Hamiltonian for the axisymmetric 
system as follows: 
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where we have followed the Lagrange multiplier subscript convention established 
earlier. Then minimizing the Hamiltonian by forming the Lagrangian of the 
Hamiltonian and evaluating the partial derivative of the Lagrangian with respect 
to the control as before we have, 
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(3.31) 


as necessary conditions for optimal control solution. Turning our attention to the 
first prolate spacecraft case (Table 3) we find that the optimal solution exhibits 
characteristics similar to the symmetric spacecraft we examined previously. The 
optimal control solution is shown in Figure 14. It displays the same bang-bang 
switching structure that was previously observed in the symmetric case. 



Figure 14 Case 1 - Axisymmetric Spacecraft Time-optimal Control 

Solution 


The quaternion and angular rate histories are shown (Figure 15 and 
Figure 16) and clearly demonstrate that the time optimal maneuver is not an 
eigenaxis maneuver. 
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Figure 15 Case I - Axisymmetric Spacecraft Time-Optimal Quaternion 

History 



Figure 16 Case I - Axisymmetric Spacecraft Time-Optimal Angular Rate 

History 

The feasibility of the solution is verified by propagating the optimal control 
solution through the state dynamics. The propagation results, shown in Figure 
17 and Figure 18, indicate that the solution does indeed drive the state from the 
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known initial conditions to the desired final conditions. The original solution 
obtained is shown in solid lines overlaid with the propagated states shown as V 
marks below. 



Figure 17 Axisymmetric Spacecraft Quaternion Solution Validation by 

Propagation 



Figure 18 Axisymmetric Spacecraft Angular Rate Solution Validation 
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The necessary conditions for optimality that we can evaluate 
directly include the switching structure obtained from the Hamiltonian 
minimization (Equation (3.31)), and the behavior of the Hamiltonian over time. 
The switching structure is shown in Figure 19 through Figure 21. The control 
solution has been overlaid to further illustrate the relationship between the 
switching function and the control solution. 



Figure 19 Axisymmetric Spacecraft x-axis Switching Structure and 

Control 
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Figure 20 Axisymmetric Spacecraft y-axis Switching Structure and 

Control 



Figure 21 Axisymmetric Spacecraft z-axis Switching Structure and 

Control 
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Inspecting the Hamiltonian for the axisymmetric case (Equation 
(3.30)) reveals no direct dependence on time. Therefore we can see that the 
Hamiltonian should be a constant over the interval under consideration. 
Additionally, forming the end-state Lagrangian as before in equation (3.24) gives 
the final value of the Hamiltonian. The evolution of the Hamiltonian over time is 
shown in Figure 22. 



Figure 22 Case 1 - Axisymmetric Hamiltonian Evolution and 

Transversality 

Our analysis of the solution indicates that it is a feasible solution to 
the time-optimal reorientation problem. Additionally, the solution meets the 
necessary conditions for optimality derived from Pontryagin’s Minimum Principle. 

Now examine the effects of increasing the length of the symmetry 
axis of a constant mass prolate spacecraft. From equations (3.29) and Table 3 
we can see that the moment of inertia about the axis cf symmetry remains 
constant while the remaining moments of inertia grow proportionally to length 
squared. This is shown in Figure 23. These increasing moments of inertia have 
two effects. First, the angular acceleration about the x & y-body axes is reduced 
proportional to the increase in moment of inertia. Physically, this is a decrease in 
the ability of a constant torque to cause an angular acceleration of an increasing 
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mass moment of inertia. However, in addition to the decreasing control 
effectiveness, we see that the switching functions of the controls about the x & y- 
axes are numerically vanishing due to the increasing moment of inertia. The y- 
axis switching function for case III illustrates this effect and is shown in Figure 25. 
Increasing the value of the control torque available is not sufficient to counter the 
vanishing switching function. 


Moment of Inertia 




Length 


Figure 23 Constant Mass Prolate Spacecraft Moment of Inertia versus 

Length of Symmetry Axis 

Therefore the time-optimal reorientation maneuver about the symmetry axis 
approaches an eigenaxis maneuver as the length of the body grows in a constant 
mass spacecraft. 
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Figure 24 Axisymmetric Reorientation Maneuvers About Symmetry Axis 

versus Symmetry Axis Length 

It is interesting to note that the structure of switching function is not 
significantly affected by the increasing moments of inertia. The switching 
function for the prolate spacecraft of case III (Table 3), where the length of the 
symmetry axis is 10 times longer than the original, displays the same switching 
structure as the original. The magnitude is however, reduced by several orders 
of magnitude. The switching function for the y-axis control of Case III is included 
(Figure 25) for comparison to that of Case I (Figure 20). 
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Figure 25 Prolate Axisymmetric Spacecraft Case III - Switching function 

b. Non-spinning Axisymmetric Reorientations of the Axis 
of Symmetry 

The case of non-spinning axisymmetric spacecraft reorientations of 
the axis of symmetry is mathematically no different from the previous section. 
The problem formulation and necessary conditions for optimality remain 
unchanged. The optimal maneuvers display the same characteristic precession 
about the eigenaxis and involve control torques about all three axes. It is 
interesting to note that the control switching structure behaves in a manner very 
similar to that observed by Bilimoria and Wie in their work on symmetric 
reorientations. 17 That is, for small angle reorientations we observe a sequential 
seven-switch structure and for large angle maneuvers we observe a sequential 
five-switch structure. The control solutions for two representative maneuvers are 
shown. In Figure 26, we see the control solution for a large angle maneuver 
about the x-axis. The representative maneuver is chosen as a 135 degree 
rotation. The five-switch-sequential structure is clearly evident. In Figure 27, the 
control solution for a representative small angle maneuver about the x-axis is 
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Figure 26 Axisymmetric Spacecraft Optimal Control Solution for X-axis 
Large Angle Rotation (135 Degree Rotation) 



Figure 27 Axisymmetric Spacecraft Optimal Control Solution for X-axis 
Small Angle Rotation (60 Degree Rotation) 
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shown. The maneuver, a 60 degree rotation results in a seven-switch-sequential 
solution. This suggests that there is a dividing point, as Bilimoria and Wie 
observed for the symmetric case, where the time-optimal solution changes from 
the seven-switch to the five-switch structure. 


c. Reorientations of the Spinning Axis of Symmetry 

Next we consider the reorientation of a spinning axisymmetric 
spacecraft. Shen and Tsiotras also considered this case where the rigid body 
was subject to only two control torques which spanned the plane perpendicular to 
the axis of symmetry. They concluded that two torques were sufficient to achieve 
a time-optimal maneuver. 18 In this section we will show that a third torque about 
the symmetry axis, if available, further reduces the objective function and is the 
true time-optimal solution. 

Consider the axisymmetric spacecraft of Figure 13, which is 
spinning about the z-axis, the axis of symmetry. Euler’ equations remain 
unchanged from equation (3.26) and are repeated here for the convenience of 
the reader: 

=/ x 03 x +(/ z -/ y )C0 y (0 z 

M 2 =/ y c6 y + (/ x -/ z )co x co z 
M 3 = / z c6 z 

Using the formulation of Shen and Tsiotras suggested by Tsiotras and Longuski 19 
the orientation of the n 3 inertial axis of the inertial frame given by 

n = (h v h 2 ,4 ) with respect to the body frame can be represented by two variables 
defined as: 

P -a 

w\ = —— w 0 =- 

1+Y 1+y 

where w i and w 2 obey the differential equations: 
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• UJ / p p \ 

W, = CO z w 2 + CO y wj/z 2 + — (1 + i/v, -W 2 J 


24 

2 

co, 


w 2 =-w z w,+ 0 ) x w,w 2 +^(l 


+ w 2 - wf 


Readers are directed to the references [18,19] for a complete derivation of this 
parameterization. Using this formulation we can state the time-optimal 
reorientation of the spin axis as follows: 


Minimize J(x{),u{),t f ) = t f -t 0 
s.t. 
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Following the numerical example of Shen and Tsiotras we establish the following 
spacecraft parameters: 


4 = 4 
4= 4 

4=2 


The Hamiltonian and necessary conditions are available in the reference and not 
repeated here. Using initial and final conditions given as: 
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x = [co 1 co 2 co 3 w 2 f 
x 0 = [0 0 -0.5 1.5 -0.5] r 
x f = [0 0 -0.5 0 Of 

we can duplicate the results of Shen and Tsiotras by setting the torque M 3 = 0. 
The results are shown in Figure 28 and Figure 29. The maneuver, defined as a 
115.38 degree reorientation of b 3 to h 3 , and is completed in 2.6142 seconds. 
Published results indicated a minimum maneuver time of 2.61 seconds. 20 



Figure 28 Angular Rate History for Constant Spin Rate Time-optimal 

Maneuver 
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Figure 29 Time Histories of Wi and W 2 for Constant Spin Time-optimal 

Maneuver 


However, if the control u 3 is available, and the boundary conditions 

are enforced such that the spin rate is allowed to vary throughout the maneuver, 
then the true time-optimal solution is found. 



Figure 30 Angular Rate History for Time-optimal Spinning Reorientation 

Maneuver 
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If spin rate of the z-axis is allowed to vary the time-optimal 
maneuver solution contains a significant change in spin rate as shown in Figure 
30. The improved time for maneuver completion is now 2.513 seconds. This is a 
reduction of 3.87% from the previously assumed time-optimal solution. The 
three-axis control time-optimal w solution is shown in Figure 31. 

The control solution obtained (Figure 32) is bang-bang in all three 
axes. As before the control solution is propagated through the state dynamics to 
verify that the solution is feasible. The original solution obtained is shown in solid 
lines overlaid with the propagated states shown as V marks below. It is clear 
(Figure 33 and Figure 34) that the spacecraft states properly propagate from the 
given initial conditions to the desired final conditions. 



Figure 31 Time Histories of w 1 and W 2 for Spinning Spacecraft Time- 

optimal Maneuver 
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Figure 32 Spinning Axisymmetric Spacecraft Time-optimal Control 

Solution 



Figure 33 Spinning Axisymmetric Spacecraft Angular Rate Solution 

Validation by Propagation 
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Figure 34 Spinning Axisymmetric Spacecraft W History Solution 

Validation by Propagation 

Minimization of the Hamiltonian with respect to the control vector 
allows us to establish the switching functions. These are shown, overlaid with 
the unity scaled control solution (Figure 35, Figure 36, and Figure 37). 
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Figure 35 Spinning Axisymmetric Spacecraft X-axis Switching Function 

and Normalized Control 

Finally, we evaluate the Hamiltonian and observe it is constant over 
time and numerically equal to -1. This is expected for the Mayer formulation of 
the cost function when the Hamiltonian has no direct time dependence. 
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Figure 36 Spinning Axisymmetric Spacecraft Y-axis Switching Function 

and Normalized Control 
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Figure 37 Spinning Axisymmetric Spacecraft Z-axis Switching Function 

and Normalized Control 



Figure 38 Spinning Axisymmetric Spacecraft Time-optimal Solution 

Hamiltonian 
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3. Numerical Considerations and Notes 

Operationally, it is often preferable to maintain a constant satellite spin 
rate about the axis of symmetry. However, we have shown that if a third torque 
is available about the axis of symmetry then the time-optimal solution will involve 
a change in the spacecraft spin rate. This change in spin rate is appreciable and 
results in a measurable time savings over the two control torque solution. 

4. Conclusions 

In this section we examined the time-optimal maneuver for an 
axisymmetric spacecraft with three independent control torques. We saw that for 
maneuvers about the axis of symmetry the relative moment of inertias have 
significant effects on the maneuvers. As the moment of inertia about the 
symmetry axis increased the maneuver approached the eigenaxis maneuver in 
both time required and spacecraft response. Additionally, we demonstrated that 
while reorientation of the spin axis is possible with two control torques spanning 
the plane perpendicular to the spin axis, the addition of a third control torque 
about the axis of symmetry further reduces the objective function and is the true 
time-optimal solution. This third control torque is in general assumed to be 
available as it was required to generate the spinning motion. 


D. ASYMMETRIC SPACECRAFT REORIENTATION MANEUVERS 

In this section we will numerically investigate the time-optimal reorientation 
of a rigid asymmetric body. Livenh and Wie 21 presented an extensive analytical 
analysis of the asymmetric reorientation problem under constant body-fixed 
torques. Additionally, the work of Proulx and Ross 22 determined an admissible 
switching structure which was elegantly illustrated by the traversal of a unit cube. 
Using this to limit the search space a combination of a genetic algorithm and 
pseudospectral method was used to obtain the optimal solution. Additionally, 
they suggested a method of evaluating the “optimality” of a solution by evaluating 
the Hamiltonian derived from the costates obtained through the Covector 
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Mapping Theorem. This method of evaluating compliance with the Minimum 
Principle is employed throughout this work. 


1. Problem Formulation 

In order to facilitate future work and extension to other control actuators 
the orientation of the asymmetric spacecraft will be represented with respect to 
an Earth centered inertial reference frame. This will require the incorporation of 
orbital velocity and reference frame transformations to properly represent the 
dynamic constraints. 

As before, for a rigid body the applied torque about the center of mass is 
equal to the time rate of change of the angular momentum. 23 


From the transportation theorem we have, 


d 

dt 





+ 


VxH 


Following the previous development we obtain Euler’s equations: 

H =4®x + (4- / y K®z 
M 2 = / y ® y + (4 - 4 )(D/o z 
M 3 = / z c6 z + (l y -4)o) x co y 


(3.32) 


The quaternion kinematics equations are unchanged from previous (equation 
(3.13)) except that emphasis is placed on the notation. 


4 = ^[ co i < 7 4 +0) 3<7 2 ] 

4 =^[«lQ 3 +^Q4-«3Ql] 

4 =^[-°V72 +“ 2 Qi +“ 3 ^ 4 ] 
4=^[-®1 C ?1-^ C ?2- 01 3^] 


(3.33) 
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As shown above the angular rates (co 1 ,co 2 ,co 3 ) are of the body with respect 

to the orbit frame. For convenience we choose to define the state of the 
spacecraft in terms of the quaternion vector and the angular rates with respect to 
the Newtonian frame. Thus the state vector is given by: 



2. Time Optimal Maneuvers 

The formal statement of the optimal control problem is as follows: 

Minimize J(x{),u{),t f ) =t f -t 0 

s.t. q = — £lq 
- 2 - 



u, <1 


The maneuver under consideration will be rest-to-rest in the orbit frame; 
however, the final angular velocity in the inertial frame will depend on the final 
attitude. This is clear from equation (3.8). The maneuver under consideration is 
an x-axis rotation typically 135 degrees. This maneuver magnitude was selected 
based on the anticipated orbital parameters of NPAST1, the potential test bed for 
later algorithms. The initial condition was chosen somewhat arbitrarily as nadir 
pointing. Then the problem initial and final conditions may be presented in 
standard form as follows: 
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x = [q i q 2 q 3 q 4 co x co y coj 

Xq=[ 0 0 0 1 0 -co 0 0] r 
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Additionally, spacecraft moment of inertias for subsequent examples have been 
selected to match the planned NPSAT1 moment of inertias and are given as 
follows: 

4 = 5 kg.m 2 
l y = 5.1 kg*m 2 
4 = 2 kg*m 2 

As before, the first step in our approach is to form and minimize the 
Hamiltonian. Recall that the Hamiltonian is a function of the state, control and 
Lagrange multipliers. 

H(X,x,u,t) = X T f{x,u) 

Since quaternion kinematics are typically written in terms of the body angular 
rates with respect to the orbit frame, a lengthy algebraic process of coordinate 
transformations is required to properly write the Hamiltonian. Once completed 
the Hamiltonian takes the following form: 


X X 

H = K< 7 4 -°V 73 + “z 9 2 +“oQ 3 ) + y K< 7 3 +C V 74 -G>z< 7 l + ) 

X X 

+ -y ( -®xQ 2 +“yQi +“ 2 Qz -®oQi ) + y-(- Q) xQi - <V h - w zQ 3 - “ 0 Q 2 ) 


(3.34) 
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In accordance with the Minimum Principle the Lagrangian of the Hamiltonian is 
formed and partial derivatives with respect to the control vector are formed which 
establish the control switching functions. 


an 

du^ 

dH_ 

du 2 

})H 

du 3 
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K 

' / 


-+ Pi - o 


■ + p 2 _ o 


^ 0 ), n 

■ —j— + \i 3 -0 


(3.35) 


The reader should note that the switching functions of the asymmetric case are 
no different than those of the axisymmetric case given in equation(3.31). This 
would correctly lead us to surmise that as the moment of inertias approach those 
previously studied the results would closely match those previously obtained. 

The control solution for the asymmetric body under consideration is shown 
in Figure 39. It clearly exhibits the structure we have come to expect. In general, 
the control is bang-bang in all three axes. There are five switches in all, 
characteristic of a large angle slew with a single-switch in the primary axis of 
rotation. 
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Figure 39 Asymmetric Spacecraft Time-optimal Maneuver Control 

Solution 

The process of solution validation begins with propagating the candidate 
solution through the state dynamics. A feasible solution must drive the 
spacecraft from its known initial state to the desired end state. The calculated 
state histories are shown figure 40 & Figure 41) with the propagation results 
(Figure 42 & Figure 43). The original solution obtained is shown in solid lines 
overlaid with the propagated states shown as V marks below. 
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Figure 40 Asymmetric Spacecraft Time-optimal Maneuver Quaternion 

History 



Figure 41 Asymmetric Spacecraft Time-optimal Maneuver Angular Rate 

History 
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Figure 42 Asymmetric Spacecraft Quaternion Solution Validation by 

Propagation. 



Figure 43 Asymmetric Spacecraft Angular Rate Solution Validation by 

Propagation. 
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The solution, confirmed feasible, is evaluated for optimality by observing 
the switching functions and Hamiltonian evolution. The switching functions, 
equations (3.35) are shown graphically with the overlaid control solution. 



Figure 44 Asymmetric Spacecraft X-axis Switching Function and Control 

Solution 



Figure 45 Asymmetric Spacecraft Y-axis Switching Function and Control 

Solution 
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Figure 46 Asymmetric Spacecraft Z-axis Switching Function and Control 

Solution. 

These figures illustrate both the Hamiltonian minimization and the system’s 
compliance with the KKT conditions. 

Inspection of the Hamiltonian reveals no direct dependence on time. Thus 
we expect a constant value Hamiltonian. The final value of the Hamiltonian is 
given by the transversality condition, equation(3.25). Taken together we expect 
a constant Hamiltonian of -1. Figure 47 shows the Hamiltonian with the 
optimality characteristics predicted. 

Finally, we compare the time-optimal solution with the eigenaxis maneuver 
once theorized as nearly time-optimal. The results, (Table 4) show a significant 
reduction in maneuver time. 
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Asymmetric Maneuver Comparison 


Eigenaxis 

6.8742 seconds 

Time-optimal 

6.1737 seconds 

Reduction 

10.19% 


Table 4 Asymmetric Maneuver Comparison 



Figure 47 Asymmetric Spacecraft Time-optimal Hamiltonian Evolution 

and Transversality 

3. Numerical Considerations and Notes 

The problem was transformed into the orbital frame as a stepping stone to 
the magnetic torque problem. Orbital position and velocity will be necessary for 
magnetic field computations. Additionally, this problem formulation incorporated 
linear scaling factors. While not strictly necessary for the problem formulation 
under consideration, incorporating scaling factors in the independent torque 
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problem allowed validation of the scaling algorithm prior to the incorporation of 
the more complicated magnetic field calculations. 

For evaluation, the torque available was reduced form 1 Newton-meter to 
0.01 Newton-meters. Without scaling, it was noted that costate estimates 
increased by 4 orders of magnitude. In this work we have found that large 
costate estimates have generally been an indication of a numerically poor 
problem formulation. Results in these cases have generally displayed poor 
accuracy and long computational times. Numerical scaling was implemented in 
the following form: 

u = k m u 
co = kjD 

T = k t t 

c/co _ cfco ( k m ' 

dt dt k, 

v 1 J 

where the overbar indicates a scaled variable. 

4. Conclusions 

In this section we have shown that the asymmetric problem displays many 
of the characteristics previously noted in symmetric and axisymmetric 
configurations. The time-optimal solution, as we might have predicted, is not an 
eigenaxis maneuver but instead is bang-bang in all control axes. The sequential 
switching structure theorized and observed by Proulx and Ross 24 was observed 
for the problem. The number of possible configurations for asymmetric 
spacecraft combined with the possible slew maneuvers limits the extent to which 
numerical analysis can be used to form general conclusions. However, the 
method employed allows us to generate the time-optimal solution to the 
asymmetric configuration, validate the feasibility of the candidate solution, and 
evaluate its compliance with the Minimum Principle. 
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IV. MAGNETIC TORQUE CONT ROL 


A. INTRODUCTION 

Magnetic torque control has been effectively used for momentum 
management in zero-momentum systems on many spacecraft systems. The 
methods are well known and flight tested. The use of magnetic torque control for 
spacecraft three-axis stabilization is less well known but the body of research is 
growing. Magnetic torque control represents a low cost method to control small 
spacecraft in reasonably low earth orbit. In this section we examine the basics of 
magnetic torque control and the solution to the time-optimal slew of spacecraft 
using magnetic torque generating devices. This problem is significantly more 
complicated then the idealized actuator problem of Chapter III. This is due to the 
resultant cross-product torque generation and varying magnetic field. It is well 
known that there are body-frame orientations where no torque can be generated 
in specific directions. 

Junkins and Turner, reference [1], discuss the magnetic time-optimal 
control of spin-stabilized spacecraft. They were able to solve the open loop 
problem for spin-axis reorientation and implement their solution on the NOVA-1 
spacecraft in 1981. 1 

B. BASIC MAGNETIC TORQUE ATTITUDE CONTROL 

The magnetic moment generated within the spacecraft, whether 
generated intentionally or inadvertently, interacts with the Earth’s magnetic field 
to produce a torque according to: 

f B = mxB 

where m is the magnetic dipole moment generated inside the spacecraft body 
and B is the Earth’s magnetic field intensity. 2 

* Magnetic dipole moment is represented as a lower case m in order to distinguish it from an 
applied torque which is commonly represented as an upper case M. 
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The Earth’s magnetic field intensity as approximated by McElvain (1962) 


is: 


B x (t) 

By(t) 

B z (t) 


cos(co 0 f)sin(/ m ) 
-cos (/ m ) 
2sin(co 0 f)sin(/ m ) 


(4.1) 


where, i m is the inclination of the satellite orbit with respect to the magnetic 
equator, R is the semi-major axis of the orbit, co 0 is the orbital angular velocity, t 
is zero at the point in the satellite orbit where the ascending node crosses the 
equator, and is the magnetic dipole strength of the Earth (circa 1975) given as: 

\i f = 7.96x1 0' 5 Wb-m 

So we see that magnetic field intensity decreases rapidly with orbital 
altitude. A typical magnetic field approximation is shown for a satellite in a 
circular orbit with inclination of 35.4 degrees at altitude of 560 Km. As expected 
the behavior in the x-z planes is harmonic (Figure 48) at the orbital frequency. 

The magnetic model used for this problem differs from equation (4.1). The 
above model is reasonably accurate as a first-order approximation but does not 
take into account the rotation of the earth. For this reason, a model with slightly 
higher fidelity was adopted. Shown below, as equation (4.2), is a model adopted 
from Wertz, Spacecraft Attitude Determination and Control. 3 This model 
assumes no orbit precession but does allow for the Earth’s rotation. 
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X' 

- 


K 

B 


Vo 


x. 

- 


cos(a) [cos(e)sin(/)-sin(e)cos(|cos(L/)]-sin(a)sin(e)sin(L/) 
-cos(e )cos (/) - sin(e )sin(/)sin(L/) 
2sin(a)[cos(£)sin(/)-sin(e)cos(|cos(L/)] + 2cos(a)sin(e)sin(L/) 


where, 


K = 7.943x10 15 -^- 
amp 

R = radius from Earth c.m. to s/c c.m. (m) 
a=co o t 

u = K> e t, co e = Earth rotational frequency (rad/sec) 
£ = 11.7° magnetic dipole tilt 
/= orbit inclination 



Figure 48 Earth Magnetic Field 


(4.2) 
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Recall Euler’s equations (repeated from previous for the convenience of 
the reader). 

=/ A+(/z-/ y K 0) z 

M 2 = / y o) y + (/, - / z K® z 
M 3 =/ z co z + ( / y 

The external torque can now be replaced by the torque generated by the 
interaction of the Earth’s magnetic field and the satellite magnetic dipole moment. 
The rotational dynamics equation of motion in vector form then becomes, 

/(6+ cox/co= mxB (4.3) 

where angular velocities are by necessity referenced to the Newtonian frame 
and, 

m 2 B z -m ? B y 

mxB = m 3 S x -m 1 S z (4.4) 

m,B - 

As written, the cross-product of equation (4.4) is meaningless since the 
dipole moments (m,) are in the body frame and the magnetic field components 

are represented in the orbit frame. Therefore, the components of the magnetic 
field which are in the orbit frame must be rotated into the spacecraft frame by: 

Be = b C°B 0 (4.5) 

where the DCM was previously defined in equation (3.5). Then the magnetic 
field in the body frame in expanded form is given by: 

B,=B Xo (<$ - $-cf 3 +ql) + 2B yo {q,q 2 + q 3 q A ) + 2 B Zo (q,q 3 + q 2 q 4 ) 

B 2 = 2 B Xq (q,q 2 -q 3 q 4 ) + B yo (4-c$-(f 3 +ql) + 2B Zo { q 2 q 3 + q 1 q 4 ) (4.6) 

B 3 = 2 B Xo (q,q 3 + q 2 q 4 ) + 2BJq 2 q 3 - q,q 4 ) + B Zq ( q 2 3 - q 2 - q\ + q\) 
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Finally, by substituting equation (4.6) into equation (4.3) we can form the 
dynamic constraint equations as follows, 


co. 


y( m 2 e 3- m 3 e 2 )+VVO z , 


«V = j{mA ~ m A ) + 


“z=f( m i e 2- m 2 e i) + ^ ro x« y - 


k Jy~K 

k '-~r 

X 

k =Lzk 

2 ly 

k 

3 L 


These represent the rotational dynamical equations of motion for our spacecraft, 
Earth magnetic field system. 


C. TIME-OPTIMAL MAGNETIC TORQUE CONTROLLED SLEW 

In this section we consider the time-optimal reorientation of a spacecraft 
with magnetic torque control. The maneuver is defined as rest-to-rest in the orbit 
frame where the initial and final states are given. 


1. Problem Formulation 

The spacecraft state is defined by its position and angular velocity. 


x ■■ 


q 

co 


The position is represented by a four-element quaternion vector, we have 
previously adopted the convention that the fourth element of the quaternion 
vector is the scalar quantity. The angular velocity is in body coordinates with 
respect to the Newtonian frame. 

We choose our control parameter as the magnetic dipole moment of the 
torque rods. Dipole moment is controlled by current flow however the response 
of dipole moment to changes in current can be considered instantaneous. Since 
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the actual torque generated is limited by the maximum dipole moment of the 
torque rod, we impose a bound on the control dipole moment. 


u = 


m 2 

m 3 


—> |/r^| < 30 A m 2 


(4.7) 


The minimum time, rest-to-rest, reorientation problem may then be stated as 
follows: 

Determine the controls [tv/, u 2 *,u 3 *]that drive the spacecraft from its initial 
rest position, given by [xq] to its final rest position given by [x f ]while minimizing 
the cost function: 

J(x(-),u(-),t f ) = t f -t 0 

where we have used the Mayer form of the cost function, subject to the following 
constraints: 

Control Constraint: The control constraint is defined in the standard form, 

h L {t)<h{u,t)<h u (t) 

then our control constraint can be written, 

-30 < u, < 30 
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Dynamic Constraints: 


Ql=^K C ?4- W 2 C ?3 +C0 3 C ?2] 

4 =^[“lQ3+^Q4-“ 3 Ql] 
4 =^[- C0 1 C ?2 +“ 2 Ql +“ 3 ^ 4 ] 
4 


m 2 B 3 ~ m 3 e 2 ) + 

X 

_N 

1 

_x 

II 

= y-(m 3 e i - m A ) + kjojtoy, 

y 

/ - / 

>< 

N 

il 

CM 

“z = y-( m A - m 2 B, ) + /c 3 co x co y , 

/- / 

k= y 

3 L 


(4.8) 


2. Solving the Optimal Control Problem 

The first step in solving the optimal control problem is to form the 
Hamiltonian. Basic format of the Hamiltonian is repeated heret, 

H(X, x , u, t) = F(x, u, t) + X T f(x, u , t) 

Since the cost functional was formulated without a Lagrange cost term the 
Hamiltonian reduces to the following. 

H{ r k,x,u,t) = X T f(x,u,t ) 

Again, since the quaternion kinematics are written in terms of angular velocity of 
the body with respect to the orbit frame, equation (4.8), a lengthy algebraic 
rotation sequence is required to write the Hamiltonian in standard form: 


t Recall that in this notation F is the Lagrange (running) cost and f are the state dynamics. 
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H = Y^x q ^ _co ^ + co z<7 2 +C0 oQ 3 ) + -| l KQ 3 +“ y q 4 -w z Q 1 + co 0 q 4 ) 

X X 

+ -^(-«xQ 2 +®yQl +“zQ 2 -“ 0 Qi) + -^(-®xQi -«yQ 2 “© zQs ~© 0 Q 2 ) + 

^ |f( W 2 e 3 - U 3 e 2 ) + *.®y®z} + V |f ( U 3 e i " ^3 ) + Mx“z 

^jy-lwi^-^ej + ^coy 

The subscripts on the Lagrange multipliers have been chosen for bookkeeping 
purposes. The control vector is defined in equation (4.7). 

Now, according to Pontryagin’s principle the control which minimizes the 
cost functional must meet the conditions we established earlier. It is however, 
important to note that not all of these conditions reveal usable information about 
the nature of the problem. 



a. Hamiltonian Minimization 

We know from previous work that a necessary condition for the 
Hamiltonian to be a minimum with respect to the control variable is that the 
partial derivative of the Lagrangian with respect to the control equals zero. In this 
case the control that satisfies this condition must also lie within the control 
constraint set. We apply Lagrange multipliers in the form, 

H(k,x,u,t, |i) = H(X,x, u, t) + \x T h{u , t) (4.10) 

Then, by inspection of equations (4.9) and (4.10), we can write the 
necessary conditions for the Hamiltonian minimization as follows. 
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h±B 2 -^B 3 +n= 0 

4 2 4 

^B 3 -^e 1+ p 2 =o (4.ii) 

x z 

h^B,-^B 2 + [i 3 =0 

4 1 L 2 


These represent useful information that can be evaluated to validate the 
candidate solution. 


b. Hamiltonian Evolution and Final Value 

The Hamiltonian evolves in accordance with the simple equation, 

H = — (4.12) 

dt 


Previously, we dealt only with Hamiltonian equations that had no specific 
dependence on time and therefore the time-rate of change was zero. In this 
case the magnetic field of the Earth introduces a time-dependence into the 

Hamiltonian. That is Therefore, the Hamiltonian is not a constant in 

the interval under consideration. 


The final value of the Hamiltonian is given by, 


..... dE ,de _ 

H[tf] + — +x >—= 0 
' dt, dt, 


where the end manifold (e) is written in the standard form, 
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e{x f ,t f ) = 


q, -sin 


Q 2 


V 2 y 


% 

q 4f -cos 


(A \ 


V 2 7 


= 0 


Cfl y +C0 o COS((|)) 

(£> Zf — co 0 sin(4>) 


Then by inspection we can see that the final value of the Hamiltonian is negative 
one. 


H[t f ]+1 =0 —> H[t f ] = -1 


Although the Hamiltonian is not a constant for the interval under consideration its 
final value represents a second numerical figure of merit to validate the optimality 
of the solution. 


3. Numerical Results 

The numerical example for this work was taken from the Naval 
Postgraduate School’s current small satellite program, “NPSAT 1.” Designed 
primarily to allow a hands-on learning experience this satellite, still in the design 
phase, will provide three-axis magnetic torque control with a pitch wheel for 
increased stabilization. The moment of inertias and orbital parameters used in 
the numerical examples were taken from NPSAT 1 preliminary designs. The 
contribution of the pitch bias wheel will be addressed in a later chapter. The 
NPSAT 1 data includes: 
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NPSAT1 Parameters 


Tab 


lx 

5 kg-m 2 

ly 

5.1 kg-m 2 

lz 

2 kg m 2 

Max Dipole Moment 

30 Amp-m 2 

Orbital Altitude 

560 km (Circular) 

Inclination 

35.4 degrees 


e 5 NPSA T1 Parameters for Numerical Simulations 


The maneuver selected for simulation is a 135 degree roll (x-axis slew). 

The time-optimal control solution for this maneuver demonstrates a 
surprisingly clean bang-bang structure. The solution, shown in Figure 49, has 10 
control switches distributed among the three axes. 



Figure 49 Time-optimal Control Solution for Magnetic Torque Problem 
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Before evaluating the optimality of the candidate solution the 
feasibility is evaluated. The control solution is propagated through a 
separate ODE 45 dynamics simulator to verify that the candidate solution 
drives the dynamic system from the initial state to the final state. The 
propagation results (Figure 50 & Figure 51) show that the control solution 
does meet the end point constraints and that the estimated states closely 
match those obtained during propagation. The original solution obtained 
is shown in solid lines overlaid with the propagated states shown as V 
marks below. 


Ol 




Tim* (ft*C) 


Figure 50 Quaternion Solution and Validation by Propagation 
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Figure 51 Angular Rate Solution and Validation by Propagation 

The quaternion and angular rate histories are shown (Figure 52 & Figure 53). 
The maneuver is clearly not an eigenaxis slew. This is evident from both the 
variation in the quaternions q 2 &g 3 and the non-zero angular rates of co 2 &co 3 . 



Figure 52 Magnetic Torque Slew Time-optimal Quaternion History 
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Figure 53 Magnetic Torque Slew Time-optimal Angular Rate History 

Next we evaluate the optimality of the feasible, candidate solution. The 
switching functions are given in equations (4.11). These are plotted overlaid with 
the scaled control solution (Figure 54, Figure 55, & Figure 56). 



Figure 54 Magnetic Torque Control X-axis Switching Function and 

Control Solution 
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Figure 55 Magnetic Torque Control Y-axis Switching Function and 

Control Solution 



Figure 56 Magnetic Torque Control Z-axis Switching Function and 

Control Solution 
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Shown are the switching functions (S,) previously defined as the partial 

derivative of the Hamiltonian with respect to the control vector. The KKT 
multiplier (Mu) is also shown plotted separately from the switching function. The 
sum of the switching function and the KKT multiplier is the definition of the 
minimization of the Hamiltonian, equation (4.11), and should be numerically 
equal to zero. Additionally, as before the switching function and control are 
related be the KKT conditions: 


U; 


maximum 
\ minimum 
singular 


S ( . <0 

S,> 0 

s,=o 


These figures clearly illustrate that the control solution meets optimality criteria 
established by the Hamiltonian minimization. 

The Hamiltonian transversality and evolution conditions are complicated 
by the varying magnetic field. The Hamiltonian and the predicted final value are 
shown (Figure 57). The numerical final value of the Hamiltonian and the 
theoretical final value differ by only 0.0646. Therefore we conclude that the 
candidate solution has met the necessary conditions for optimality. 



Figure 57 Magnetic Torque Solution Hamiltonian and Final Value 
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Additionally, to validate our hypothesis that the time-dependence of the 
Hamiltonian was due to the varying magnetic field the problem was evaluated 
with a constant magnetic field approximation. The Hamiltonian for the case of a 
constant magnetic field is shown in Figure 58. The Hamiltonian for this case has 
lost its apparent dependence on time and settled to a value that is numerically 
close to the value of -1 that was predicted. 


0 

-0.5 

-1 

-1.5 

-2 

o c 

Constant B-Field Hamiltonian 
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100 200 300 400 500 


Figure 58 Constant Magnetic Field Approximation Hamiltonian Solution 
4. Numerical Considerations and Scaling 

The linear scaling used throughout these algorithms was introduced in the 
previous section. In this section scaling was also added for the Earth’s magnetic 
field in the form: 

B, = k B B, 

The goal of the scaling is to bring all numerical values seen by the optimization 
solver into the same order of magnitude. Scaling values were adjusted from an 
unsealed solution to improve the quality of the solution and then readjusted as 
necessary. Proper scaling reduced computation time and improved the accuracy 
of the solution. Additionally, the switching functions of properly scaled problems 
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behaved more closely to theoretical predictions indicating that results were 
improved by proper scaling. 


Scaling Effects 

Maneuver Time (before 

scaling) 

230.0845 seconds 

Maneuver Time (after 

scaling) 

271.1564 seconds 

Error 

15.15% 


Table 6 Effects of Scaling on Solution Fidelity 


5. Conclusions 

In this section the open loop time-optimal control for a magnetic torque 
controlled asymmetric spacecraft was determined. The candidate solution was 
determined by propagation to be a feasible solution to the problem. The 
optimality of the solution was validated through an analysis of the Hamiltonian 
minimization, switching functions and the behavior of the Hamiltonian. The 
hitherto unseen variation of the Hamiltonian over time was theorized to be 
caused by the time dependence of the Earth’s magnetic field. When this 
dependence was eliminated the Hamiltonian returned to the constant values we 
have seen previously. Therefore, we conclude that the solution is feasible and 
meets the necessary conditions for optimality. 
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V. NPSAT 1 CONTROL SYSTEM 


A. INTRODUCTION 

NPSAT 1 is a small satellite design project currently in work at the U.S. 
Naval Postgraduate School (Figure 59). It was conceived as a three-axis 
stabilized magnetic torque controlled satellite. Later in the design process, a 
pitch bias wheel was added to improve stability and reliability. In this section we 
explore the time-optimal reorientation of this small asymmetric satellite where 
both the three torque rods and the pitch wheel are available as torque generating 
devices. 



Figure 59 NPSA T1 Conceptual Image courtesy of Dan Sakoda 
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B. PROBLEM FORMULATION 


In this case the spacecraft state is defined by its position, angular velocity 
and the angular momentum of the pitch wheel. 


x = 


Q 

co 

h 


The position is represented by a four-element quaternion vector; we have 
previously adopted the convention that the fourth element of the quaternion 
vector is the scalar quantity. The angular velocity is in body coordinates with 
respect to the Newtonian frame. The pitch wheel is assumed to be aligned with 
the spacecraft #2 principal axis. This assumption simplifies the formulation of the 
rotational dynamics equations. The orbital parameters are as given in Table 5. 


The quaternion kinematics equations are well known and shown in 
equation (3.33). These equations are unchanged by the addition of the pitch bias 
wheel. The rotational dynamics however, are now the result of the external 
torque generated by the interaction of the magnetic field and the torque rods in 
addition to the pitch wheel angular momentum and torque. 

The time rate of change of angular momentum can be expressed as 1 , 


M ext = 


d n 


d n 

—H q 

— 

—H q 

[dt s \ 

N 

[dt s \ 


+ (ox H, 


(5.1) 


where H s is the total angular momentum of the spacecraft-wheel system and is 

expressed in the body frame. Then by assuming that the wheel’s center of mass 
is collocated with the spacecraft center of mass we can express the system 
angular momentum as, 

H = l B % b + l w N 65 w (5.2) 


where l B &l w are the body and wheel moments of inertia respectively, and the 
angular velocity of the wheel with respect to the Newtonian frame is given by, 



Then the angular momentum of the system can be written as, 


By defining, 


H q 


l B N (D b 


+ i w N <» b 


+l w b ti w 


J + l]/V 


we have, 

H s = J N (5 b + l w b 65 w 


(5.3) 


This result is easily derived under the assumptions stated. Kane 2 provides a 
detailed derivation to show that this result holds for any configuration. 

Then, referring to equation (5.1) we have, 

M ext = J N (b b + l w b cb w + N (5 b x(J N <5 b +l w b (5 w ) 


By defining, 



•w 


b ti w 



we can write, 

M ext = J N (5 b + h w + N (5 b xJ N (5 b + h w 


If we allow M ext to be the sum of disturbance torques and the torque generated 
by the interaction of the Earth’s magnetic field and the magnetic torque rods, 
setting the disturbance torque to zero we can write, 


Jco +(3x Jco = mx B-h w ~(3xh w 
where, 



1 

' O 

_i 


1 

o 

1 _ 

h w - 

h w 

> "w — 

h w 


l 

O 

i _ 


1 

o 

_ 1 


and co = N ti b 


Carrying out the cross products and rotational transformations previously defined 
gives the following result: 
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JiCO x - ( J 2 - J 3 K«z = m 2 e 3 - m 3 B 2 + CO A 
J 2 co y - (J 3 -J, )co/o z = m 3 B ; - m,B 3 - h w 
J 3 co z -(J, - J 2 )co/o y = m,B 2 - m 2 B 1 -co A 

where the Earth’s Magnetic Field is in the spacecraft body frame and 

was previously defined by equation (4.6). The four-element control vector is 
chosen as: 

A" 

m 7 

u= 2 (5.4) 

m 3 

_u w _ 

Box constraints are imposed on the control elements based on the physical 
limitations of the control elements selected. The components under 
consideration in this model have the physical characteristics shown (Table 7). 


Component Characteristics 

Torque Rods 


Maximum Magnetic 

Dipole Moment 


30 Amp*m 2 

Pitch Wheel 


Maximum Angular 

Momentum 


18 N*m*sec 

Pitch Wheel 


Maximum Torque 


30 mN*m 


Table 7 NPSAT1 Simulation Component Characteristics 


Pitch wheel angular momentum limits are imposed as a state constraint. Other 
values are imposed as control constraints. 

C. TIME-OPTIMAL NPSAT 1 SLEW 

In this section we consider the time-optimal reorientation of an asymmetric 
spacecraft controlled by a combination of magnetic torque rods and a pitch 
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momentum wheel. The maneuver is defined as rest-to-rest in the orbit frame 
where the initial attitude and attitude rates are known. The final attitude and 
attitude rates are determined by the eigenaxis of rotation and the rotation angle. 
The angular momentum of the wheel is the final state variable. For this variable 
and control combination there are four possible options. First, the initial wheel 
angular momentum may be left unspecified as an optimization variable to be 
determined. Then, if wheel torque does not equal zero, the algorithm selects the 
initial wheel speed and control history for the minimum time maneuver. This is 
the case that is numerically evaluated below. The wheel speed can also be fixed 
at the end points. In this case a non-zero toque limit will achieve the time optimal 
solution within the constraints provided. Under the condition of zero torque, 
wheel speed free, the algorithm will determine the optimal, constant wheel speed 
for the desired maneuver. Finally, if wheel speed and torque are set to zero the 
results match those previously obtained for the magnetic torque control section. 
It is important to note that these changes in boundary conditions do not affect our 
ability to minimize the objective function and obtain valid solutions. 

1. Problem Statement 

Determine the controls [u*,u 2 *,u z *,u A *] that drive the spacecraft from its 
initial rest position, given by [xq] to its final rest position given by [x f ]while 
minimizing the cost function: 


where we have used the Mayer form of the cost function, subject to the following 
constraints: 

Control / State Constraint: The control and state constraints are defined in 
the standard form, 

h L {t)<h{u,t)<h u {t) 
then our control constraints are written, 
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-30 < u t <30 (Amp-m 2 ) / = 1,2,3 
-30x10 3 < l/ 4 < 30x10“ 3 (N m) 

and the state constraint is given by: 

-0.18 <x 8 <0.18 (N-m sec) 


Dynamic Constraints: 


9l =^KQ4-^Q 3 + 0) 3 C ?2] 

Qz = - co 3 Qi ] 

+°V7i +© 3 <7 4 ] 

“x = -7- (m 2 B 3 - m 3 B 2 -a z h w ) + kp y a z , 

J x 

“y = ")"( - m v B 3 + K ) + M#,. 

J y 

®z = J-(^e 2 - m 2 B ; + G>A)+MxC0y. 




J y~ J Z 

J, 


k _J Z ~J X 

2 ~ 

, 


(5.5) 




h w = u * 
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2. Solving the Optimal Control Problem 

The Hamiltonian is given by, 


X X 

H = y (® x q 4 -a y q 3 +co zQ 2 +“oQ 3 ) + y (co x g 3 +co y q 4 -to z q A +co 0 q 4 ) + 
X X 

y(-®xQ 2 +“/7i +to 2 q' z -to 0 qO+yHv7i -co y g 2 -co z g 3 -co 0 g 2 ) + 


1 


1 


-T ( U 2 e 3 - U A -®z^ ) + ¥,®z + V T ( U 3 S 1 - "A + '“J + + 


J 




V j — (<"A - s, + CO A) + k 3 co x (o y [ + \ u 4 


J. 


The subscripts on the Lagrange multipliers have been chosen for bookkeeping 
purposes. The control vector is defined in equation (5.4). 


Following the solution method previously established we begin by 
minimizing the Hamiltonian subject to the control constraint set. Then the 
necessary conditions for the Hamiltonian minimization are as follows: 




X 


-^b-^b 3 + M =0 


J 


J 


J x J z 

^e,-^e 2 +n 3 = o 

y u x 


K 

j. 


~ + Xh + M-4 - 0 


muj 


(5.6) 


These relationships will be evaluated to validate the candidate solution. 

We saw that in the case of the magnetic torque problem the magnetic field 
of the Earth introduced a time-dependence into the Hamiltonian. Therefore, the 
Hamiltonian was not a constant in the interval under consideration. For the 
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current problem the Earth’s magnetic field is still a factor in the Hamiltonian. 
However, in this case, we will see that the effect of the magnetic field is reduced 
by the introduction of the pitch wheel, and the time rate of change of the 
Hamiltonian is reduced significantly. 

The final value of the Hamiltonian is again given by, 


dE t de n 
H[t f ] + —+W — =° 


at 


dt f 


where the end manifold (e) is written in the standard form previously established, 

e{x f ,t f )= 0 

Then by inspection we can see that the final value of the Hamiltonian is negative 
one. 

H[f,]+1=0-»H[f,] = -1 


3. Numerical Results 

The numerical example for this work was again taken from the Naval 
Postgraduate School’s current small satellite program, “NPSAT 1.” The moment 
of inertias and orbital parameters were previously established in Table 5 and are 
assumed to be the system moments of inertia. The control parameters are given 
in Table 7. 

The maneuver selected for simulation is a 135 degree roll (x-axis slew) 
and is a rest-to-rest maneuver in the orbit frame. The initial and final pitch wheel 
angular momentums are free within the state bounds as an optimization 
parameter. The time-optimal control solution is shown in Figure 60 and Figure 
61. 
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Figure 60 NPSAT1 Torque Rod & Wheel Time-optimal Control Solution 

(V 



Figure 61 NPSAT 1 Pitch Wheel Torque Time-optimal Control Solution (2) 


There are several interesting characteristics to this solution not the least of 
which is an apparent singular arc in the switching structure for the pitch wheel 
torque ( u 4 ). Before evaluating the optimality of the candidate solution the 

feasibility is evaluated. The control solution is again propagated through a 
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separate ODE 45 dynamics simulator to verify that the candidate solution drives 
the dynamic system from the initial state to the final state. Previous control 
solutions were propagated using a linear interpolation. In this case a piecewise 
cubic hermite interpolating polynomial (pchip) produced more accurate results. 
The propagation results figure 62, Figure 63, and Figure 64) show that the 
control solution does meet the end point constraints and that the estimated states 
closely match those obtained during propagation. The original solution obtained 
is shown in solid lines overlaid with the propagated states shown as V marks 
below. Therefore we conclude that the control solution is feasible. 

The state histories are shown figure 65, Figure 66, and Figure 67). It 
comes as no surprise that the time-optimal maneuver is not an eigenaxis 
maneuver. This is evident from both the variation in the quaternions q 2 &g 3 and 
the non-zero angular rates of co 2 &co 3 . 
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Figure 63 Angular Rate Solution and Validation by Propagation 



Figure 64 Pitch Wheel Momentum Solution and Validation by 

Propagation 
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Figure 66 


NPSA T1 Slew Time-optimal Angular Rate History 
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Figure 67 NPSAT1 Slew Time-optimal Pitch Wheel Angular Momentum 

History 

Next we evaluate the optimality of the feasible, candidate solution. The 
switching functions are given in equations (4.11). These are plotted overlaid with 
the scaled control solution (Figure 68, Figure 69, Figure 70, and Figure 71). 



Figure 68 NPSAT 1 Control Switching Function and Control Solution (1) 
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Figure 69 NPSAT1 Control Switching Function and Control Solution (2) 



Figure 70 NPSAT 1 Control Switching Function and Control Solution (3) 
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Figure 71 NPSAT1 Control Switching Function and Control Solution (4) 


Shown are the switching functions (S,) previously defined as the partial 

derivative of the Hamiltonian with respect to the control vector. The KKT 
multiplier (Mu) is also shown plotted separately from the switching function. The 
sum of the switching function and the KKT multiplier is the definition of the 
minimization of the Hamiltonian, equation (4.11), and should be numerically 
equal to zero. As before the switching function and control are related by the 
KKT conditions: 


u. 


maximum 
■ minimum 
singular 


S ( . <0 

Si >0 

Sj -0 


The Hamiltonian for this problem is shown in Figure 72. As we alluded to 
earlier the time dependence observed in the magnetic torque control problem 
appears to have been eliminated. In fact, it has been reduced to the extent that it 
is no longer visible. By constraining the pitch wheel to zero torque and zero 
angular momentum the previous magnetic torque results are obtained. 
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Hamiltonian Evolution 



Figure 72 NPSAT1 Torque Rod & Wheel Problem Hamiltonian Evolution 

These figures clearly illustrate that the control solution meets optimality 
criteria established by the Hamiltonian minimization. A closer look at the 
Hamiltonian minimization with respect to the pitch wheel torque is given in Figure 
73. The switching function appears singular during the period in which the pitch 
wheel torque is zero. 



Figure 73 NPSAT 1 Control Switching Function for Wheel Torque (Close 

up) 
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4. Numerical Considerations and Scaling 

We continue with the linear scaling established previously. In this section 
scaling was also added for the pitch wheel parameters in the form: 

h=k h h 

h =^h 

k t 

As before the goal of the scaling is to bring all numerical values seen by the 
optimization solver into the same order of magnitude. Scaling values were 
adjusted from an unsealed solution to improve the quality of the solution and then 
readjusted as necessary. 

5. Conclusions 

In this section the open loop time-optimal control for an asymmetric 
spacecraft equipped with three magnetic torque rods and a pitch wheel was 
determined. The candidate solution was determined by propagation to be a 
feasible solution to the problem. The optimality of the solution was validated 
through an analysis of the Hamiltonian minimization, switching functions and the 
behavior of the Hamiltonian. Therefore, we conclude that the solution is feasible 
and meets the necessary conditions for optimality. 

By allowing the use of the pitch wheel as a spacecraft control, the time 
required for the reorientation maneuver was significantly reduced. In the original 
configuration, with only magnetic torque control available, the time required for 
the optimal maneuver was 271.1564 seconds. In the new configuration the same 
maneuver required only 141.9933 seconds. This represents a 47.6 percent 
reduction in time required for the same maneuver. 
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VI. CONTROL MOMENT GYRO CONTROL SYSTEMS 


A. INTRODUCTION 

Control moment gyros (CMGs) are well known for their large torque 
generating capability. However, in the past, because of their mechanical 
complexity, cost and weight they have been restricted to large satellites with high 
agility requirements. Recent research in the design of smaller, less expensive 
CMGs has created a renewed interest in CMG spacecraft control. However, the 
primary focus of this large body of work has been singularity avoidance 12 3 4 . In 
this section we calculate the time-optimal solution for an asymmetric spacecraft 
using the much studied four-CMG pyramid configuration. In addition to 
completing a time-optimal maneuver we will show that the problem formulation 
generates a singularity free solution. 

B. CONTROL MOMENT GYRO BASICS 

A control moment gyro contains a flywheel which spins at a constant rate. 
The spin axis of the flywheel is connected to a gimbal which allows reorientation 
of the spin axis. Therefore, the direction of the angular momentum vector can be 
changed with respect to the spacecraft body frame by gimballing (Figure 74). 


Single-gimbal Control Moment Gryo 


Torque Vector 


Gimbal M 





Gyro Motor 


Flywheel 


Angular Momentum Vector 


Figure 74 Single-Gimbal Control Moment Gyro (After Ref.[5]) 
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For single-gimbal control moment gyros (SGCMG), the spin axis reorientation is 
restricted to a plane perpendicular to the gimbal axis. The advantage of CMGs is 
that small gimballing results in a large control torque on the spacecraft. This 
resultant torque is orthogonal to both the gimbal and spin axes. 6 The torque from 
the CMG is given by: 

^cmg — hx§ ( 6 . 1 ) 

where h is the angular momentum vector and has units of Newton-meter- 

seconds, 5* is the gimbal angle rate with units of radians per second. This 
torque amplification property makes CMGs desirable for applications that require 
spacecraft agility. 


C. PROBLEM FORMULATION 

Following a similar development to that used in Chapter V, recall that the 
time rate of change of angular momentum can be expressed as, 


M ext = 


d r , 


d 

— He 

— 

— He 

[ dt s \ 

N 

[dt s \ 


+ (5 X Hr 


( 6 . 2 ) 


where H s is the total angular momentum of the spacecrat-CMG system and is 

expressed in the spacecraft body frame. Assuming that the CMG center of mass 
is collocated with the spacecraft center of mass we can express the system 
angular momentum as: 

n - l B (O + i CMG co 


where l B and l CMG are the spacecraft and CMG system moment of inertias 

respectively and the angular velocity of the CMG with respect to the Newtonian 
frame is given by, 

n (5 cmg = w 63 b + b ti CMG 

Then the angular momentum of the system can be written as, 
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Once again, by defining, 


J — Ib + ICMG 


we have, 

H s = J% b + l CMG b (5 CMG 

Then referring to equation (6.2) we have, 

M ext = J N d5 b + l CMG b <b CMG + % b x(J N & b + l CMG b (3 CMG ) 

Finally, by defining, 

— I b (f) G ^ G ~ h 

~ CMG UJ _ ' CMG 

we obtain the following: 

M ext =J% b + h CMG+ % b x J% b + h CMG 

If we allow w co fa = co and M ext to be all disturbance torque, then setting 
disturbance torque equal to zero we write, 

J(5 + co x J(5 = -h CMG - (3 x h CMG (6.3) 

The CMG angular momentum vector (/?)* is a function of the gimbal 
angles (5 ) in the body frame and for multiple CMG configurations, a function of 
the configuration. Here we consider the four-CMG pyramid configuration (Figure 
75). Each face of the pyramid is inclined from the horizontal by a skew angle 
(P ). The four CMGs have gimbal axes orthogonal to the pyramid faces and so 

are constrained to gimbal on the faces of the pyramid. We have selected a skew 
angle, p = 54.73 degrees. For CMGs with equal angular momentum about the 


d_r 

dt ''CMG 


* We have adopted the notation lower case ‘h’ as CMG angular momentum and subsequently 
dropped the subscript. 



spin axis this configuration results in a nearly spherical momentum envelope. 
Additionally, this configuration has been extensively studied in the literature. 7 


Reference 

Frame 


CMG #2 
Gimbal Axis 

S, 

Pitch Axis 

CMG #1 
Gimbal 
Axis 


Figure 75 Pyramid Mounting Arrangement of Four Single Gimbal CMGs 

(from Ref.[8]) 

Then referencing Figure 75, we can write the total CMG angular 
momentum expressed in the body frame as: 

-cpsinSJ [ -cosS 2 1 [cpsinSgl cosS 4 

h = cos8 1 + -c(3sin8 2 + -cosS 3 + cpsinS 4 (6.4) 

spsinS^ sp sinS 2 J |_sp sin5 3 J |_sp sin5 4 

where cp = cosp and sf> = sin p . The angular momentum magnitude is set to 
unity without loss of generality {h 0 = 1). Then the time of rate of change of 
angular momentum may be written as: 

-cpcosS 1 sinS 2 cpcosS 3 

-sinS 1 -cp cosS 2 sin8 3 

sp cosS! spcos8 2 spcos8 3 

A 
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The matrix A is defined from equation (6.4) as, 


and is a 3 xn Jacobian matrix where n is the number of CMGs in the 
configuration. Expanding equation (6.3) showst, 

~ K - ^4K<»z = -g>24 +“ 3 4 

J y“y - ( J z - J >x®z = ~K + ®A 

~( J x~ J y )°¥°y = “4 - 4 + ®2 ^ 

where the values of fyare established in equation (6.4) and the values of h : are 

established in equation (6.5). A lengthy direct substitution allows us to write the 
time rate of change of angular velocity in the standard form: 

® i = f{x i ,u i ) 


where the state vector and control vector are defined as: 

state: [q,co,6] r e M 11 
control: [u] r e K 4 

D. CMG TIME-OPTIMAL SLEW MANEUVERS 

In this section we consider the time-optimal reorientation of an asymmetric 
spacecraft controlled by control moment gyros. In order to simplify the 
formulation the spacecraft is assumed to be in inertial space. The maneuver is 
defined as rest-to-rest in the inertial frame where the initial attitude and attitude 
rates are known. The final attitude and attitude rates are specified. The gimbal 
angle of the CMG, the final state variable is left free as an optimization variable. 


t This derivation is also found in reference [6] where it is used in the development of a 
feedback control law. 
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1. Problem Statement 

Determine the controls [u^*,u 2 *,u 3 *,u 4 *] that drive the spacecraft from its 
initial rest position, given by [xq] to its final rest position given by [x f ]while 
minimizing the cost function: 

J(x(-),u(-),t f ) = t f -t 0 

where we have used the Mayer form of the cost function, subject to the following 
constraints: 

Control / State Constraint: The control and state constraints are defined in 
the standard form, 

h L {t)<h{u,t)<h u (t) 
then our control constraints are written, 

-1 <U;< 1 rad/sec / = 1,2,3,4 

and the state constraint is given by: 

-7i < 8, < 71 rad / = 1,2,3,4 

Physically, these equate to a CMG that is capable of full 360 degree rotation 
about its gimbal axis at a rate of 1 radian per second. With a unit angular 
momentum then from equation (6.1) the maximum torque output from the CMG is 
1 Newton-meter. 
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Dynamic Constraints: 


q = —£lq 
- 2 =- 

co x =^—co y co 2 --^ 

. J y — Jy hn 

co y = ^-^co x to 2 


_h. 

Jy ' X " y Jy 


CO, = ———CO „co„ — — 


7 K^-coA) 

J x 

—(C0 z /7i-C 0 x /? 3 ) 
J y 

— (cO x /? 2 -CO y /7i) 


( 6 . 6 ) 


8/ = U; /= 1.4 

where the quaternion kinematics equations shown have been previously defined. 
The quantity /?, is defined in equation (6.5) in terms of the state 5 and the control 

and the quantity /?, is defined in equation (6.4) in terms of the state 5 . 


2. Solving the Optimal Control Problem 

As before, the first step in solving the optimal control problem is to form 
the Hamiltonian. Basic format of the Hamiltonian is repeated here*, 

H(k, x, u, t) = F(x, u, t) + X T f(x, u, t) 

Since the cost functional was formulated without a Lagrange cost term the 
Hamiltonian reduces to the following. 

H{X,x,u,t) = X T f{x,u,t) (6.7) 

Substituting equation (6.6) into equation (6.7) gives the Hamiltonian for the CMG 
spacecraft system. 


i Recall that in this notation F is the Lagrange (running) cost and f are the state dynamics. 
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" = W 


-aq 
o =— 

v. y 

y 


+ X 


J y~ J Z 


V 




Vz-j-jK ^-“A) 


+ 


J y J y J y 

f - Jx . Jy Q)/O y -4-^M ~<V*i) 


V 


V J . 


4 4 


+ 


+ 


( 6 . 8 ) 


y 


V" 


Recall from equation (6.7) that the Hamiltonian is written in terms of the 
Lagrange multipliers, state vector, control vector and time. In order for equation 
(6.8) to be rigorously correct substitutions from equation (6.4) and equation (6.5) 
are required. However, this lengthy substitution is omitted. 

Since the control is constrained, Hamiltonian minimization is accomplished 
by adjoining the control constraint equations to the Hamiltonian as, 

H(k,x,u,t,\i) = H(k,x, q t) + \i T h(u, t) (6.9) 


Then, by differentiation of equation (6.8), with respect to the control vector and 
substituting from equation (6.4) and (6.5), we can write the necessary conditions 
for the Hamiltonian minimization as follows. 


— =—cos(3 cosS, + -^sin8 1 sin (3 cosS^A,^ +p 8i =0 
dU 1 J x Jy J z 


dH 

du 2 

dH 

du^ 


k sin52 + ^cos p C0S 5 2 _ k sin p C osS 2 + ^ + - 0 


J, 


-cosp COSSg 


J„ 


■sinS q 


K 

j. 


■sin(3 cos8 3 + \ +p =0 


-sin8 4 —^cos(3cos 8 4 sin (3cos8 4 + \ +(i 5 =0 


du. J „ 


J,, 


J' 


MU: 


( 6 . 10 ) 


These conditions will be evaluated to validate the candidate solution. 
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Inspection of the Hamiltonian equation reveals no direct dependence on 
time. Therefore the time rate of change of the Hamiltonian is zero and we expect 
a constant valued Hamiltonian for the interval under consideration. 

The final value of the Hamiltonian is given by, 

dE t de 

H[t f ] + —+t> — = 0 
f dt f dt f 

where the end manifold ( e) is written in the standard form previously established, 

e{x f ,t f )=0 

Then by inspection we can see that the final value of the Hamiltonian is negative 
one. 


H[t f ]+1 =0 —> H[t f ] = -1 (6.11) 

Thus we expect a Hamiltonian with a constant value of -1 over the interval of the 
maneuver. 


3. Numerical Results 

For this numerical example we used the asymmetric moment of inertia of 
the NPSAT 1 small satellite design previously established (Table 5) and 
assumed to be the system moments of inertia. The orbital parameters were 
neglected since the model was assumed in inertial space. The physical 
parameters for the CMG are summarized in Table 8. 


NPSAT 1 SIMULATED CMG SYSTEM 

Gimbal Rotation Range 

360 deg 

Maximum Gimbal Rotation Rate 

1 rad/sec 

Maximum Angular Momentum 

1 N-m-sec 

Maximum Output Torque 

1 N-m 


Table 8 Simulated CMG Characteristics 
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The maneuver selected for simulation is a 135 degree roll (x-axis slew) 
and is a rest-to-rest maneuver in the inertial frame. The initial and final CMG 
gimbal angles are free within the state bounds. Therefore, we can write the end 
point conditions as follows: 


x, 


x = [q,,q 2 ,q 3 ,<Q 4 ^ § 3 >< §4 I r 

x 0 = [0,0,0,1,0,0,0,8^ ,8 % \ ,8 4o f 


sin 


135 


,0,0,cos 


135 


-i7 


,0,0,0,8 1 ,8 2 , 8 3 ,8 4 


The time-optimal control solution is shown in Figure 76. Our engineering 
intuition would lead us to expect a bang-bang solution for a time-optimal 
maneuver. The control solution indicates that the solution is a bang-bang 
response from all four CMGs with 5 control switches. This is a surprisingly clean 
solution considering the complexity of the CMG dynamics derived above. 

As we have seen previously, the feasibility of the candidate solution is 
evaluated propagating the control solution through a separate ODE 45 dynamics 
simulator to verify that the candidate solution drives the dynamic system from the 
initial state to the final state. A linear interpolation was used in the propagation 
sub-routine. The propagation results (Figure 77, Figure 78, and Figure 79) show 
that the control solution does meet the end point constraints and that the 
estimated states closely match those obtained during propagation. The original 
solution obtained is shown in solid lines overlaid with the propagated states 
shown as V marks below. Therefore we conclude that the control solution is 
feasible. 
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CMG Control History 



Time (sec) 


Figure 76 CMG Time-optimal Control Solution 



Figure 77 CMG Solution Quaternion History Validation by Propagation 
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Figure 78 CMG Solution Angular Rate History Validation by Propagation 



Figure 79 CMG Solution Gimbal Angle History Validation by Propagation 
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The state histories are shown (Figure 80, Figure 81, and Figure 82). Once 
again we see that the time-optimal maneuver is not an eigenaxis maneuver. This 
is evident from both the variation in the quaternions g 2 &g 3 and the non-zero 

angular rates of a> 2 &co 3 . Additionally, the gimbal angle history demonstrates 



Figure 81 CMG Slew Time-optimal Angular Rate History 
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the state is reasonably well behaved and within the state constraints [-71,71] 

imposed in the problem formulation. The time required to complete the 
maneuver is 4.24 seconds. This, as expected, is faster than the previous 
idealized actuator solution for this moment of inertia configuration (see Chapter 
III). The torque available from the four actuators is higher than the previous 
idealized actuator toque numerical example and therefore we expect faster 
maneuvering. So our engineering judgment leads us to believe that the solution 
is correct. 

Next we evaluate the optimality of the feasible, candidate solution. The 
switching functions are given in equations (6.10). These are plotted overlaid with 
the control solution (Figure 83 through Figure 86). Shown are the switching 
functions (S, ) previously defined as the partial derivative of the Hamiltonian with 

respect to the control vector. The KKT multiplier (Mu) is also shown plotted 
separately from the switching function. The sum of the switching function and 
the KKT multiplier is the definition of the minimization of the Hamiltonian, 
equation(4.10), and should be numerically equal to zero. 


120 














Figure 83 CMG Control Switching Function and Control Solution 



Figure 84 CMG Control Switching Function and Control Solution (2) 
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Figure 85 CMG Control Switching Function and Control Solution (3) 



Figure 86 CMG Control Switching Function and Control Solution (4) 
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As before the switching function and control are related by the KKT conditions: 


u. 


maximum 
• minimum 
singular 


§ <0 
S,> 0 
Sj -0 


These figures clearly illustrate that the control solution meets the necessary 
conditions for optimality established by the Hamiltonian minimization. 

The Hamiltonian transversality and evolution conditions were defined in 
equations (4.12) and (6.11). The computed Hamiltonian is shown in Figure 87. It 
is clear that the Hamiltonian is numerically well behaved and meets the 
necessary conditions for optimality. 



Figure 87 CMG Time-optimal Slew Solution Hamiltonian 


Based on an analysis of the necessary conditions and engineering 
judgment the solution obtained appears optimal for the system under 
consideration. 
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4. Numerical Considerations 

The principal problem associated with CMG systems for spacecraft 
attitude control is the well known geometric singularity conditions This refers to a 
condition in which no torque is generated for a commanded control torque in a 
certain “singular” direction. Mathematically, this occurs when the A matrix 
defined in equation (6.5) is singular or rank deficient. By solving for the optimal 
control vector u* to minimize an objective function we have required that the 

control solution satisfy equation (6.5) as part of the constraints in the problem 
formulation. As a result the possibility of the matrix A being singular is 
eliminated. The condition number of A is an indication of how close this matrix is 
to being singular^. Condition number is defined as the ratio of the largest to the 
smallest singular value of the matrix. The condition number of an identity matrix 
is one and the condition number of a matrix approaches infinity as the matrix 
becomes singular. The condition number of the A matrix is shown in Figure 88. 
It is clear that the matrix is well behaved throughout the maneuver. 



Figure 88 Condition Number of CMG Control Solution Jacobian Matrix 
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The far-reaching effects of this singularity free solution will be discussed in the 
next chapter. 


5. Conclusions 

In this section the open loop time-optimal control for an asymmetric 
spacecraft equipped with four CMGs arranged in a pyramid configuration was 
determined. The candidate solution was determined by propagation to be a 
feasible solution to the problem. The optimality of the solution was validated 
through an analysis of the Hamiltonian minimization, switching functions and the 
behavior of the Hamiltonian. Therefore, we conclude that the solution is feasible 
and meets the necessary conditions for optimality. Additionally, it was shown 
that the control solution is completely free of singularities. This leads us to 
conclude that if the control solution can be computed in real-time, then singularity 
avoidance in control moment gyro actuator based systems is unnecessary. 
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VII. CLOSED LOOP CONTROL 


A. INTRODUCTION 

Classical closed-loop control techniques are, in general, ill suited for slew 
maneuver time-optimal control. This is due in part to the nonlinear nature of the 
problem and the inability of classical control techniques to generate time-optimal 
solutions. Additionally, control solutions must be resolved into actuator input 
parameters to generate the desired torque. This method, which requires 
inverting the actuator dynamics equation, has led to singularity problems that 
continue to plague engineers as evidenced by ongoing researches,4, 5 . 

This chapter addresses the computation of real-time optimal controls for 
spacecraft slew maneuvers. Prior research efforts were focused on minimizing 
deviations from a nominal optimal trajectory. These neighboring optimal control 
(NOC) laws, presented by Bryson and others 6 ’ 7 have been used in a variety of 
applications. 

More recently a two degree-of-freedom, non-linear optimal control system 
architecture was suggested by Strizzi, Yan, et al. 8 ’ 9 This concept relied on the 
computational power of pseudospectral methods to develop optimal controls and 
NOC laws. 

In this work optimal control solutions to the non-linear plant are generated 
to implement a sampled-data feedback control law.i° Previous optimal solutions 
are used as guess for subsequent re-optimized trajectories and significant 
reductions in computation time are demonstrated. 

B. CLASSICAL CLOSED LOOP CONTROL 

A basic diagram of a classical closed loop spacecraft attitude control 
system is shown in Figure 89. In this type of model sensors determine the 
attitude and attitude rates. These are in some way compared with desired 
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Motion 




Figure 89 Classic Closed Loop Spacecraft Attitude Control System 

attitude and rates and the difference or error signal is used to determine a 
desired control torque to reduce the error signal to zero. Typical approaches for 
control command generation involve using Euler angle errors for small rotations 
and either a direction cosine error matrix or a quaternion error vector for large 
angular rotations. 11 In either case proportional plus derivative gains are imposed 
in order to achieve desired performance. This type of command generation has 
two inherent problems. First, the command generation performs no better than 
an eigenaxis maneuver and in fact will seldom perform so well. Second, 
generating a desired control torque for use as an actuator input has led to 
numerous singularity problems of which the best known is the control moment 
gyro problem. 


1. Error Signal Command Generation 

In previous sections we have seen that the time-optimal solution is a 
function of the magnitude of the maneuver and the spacecraft moment of inertia. 
While error signal command generation is capable of generating the eigenaxis 
path the command signal is not the familiar bang-bang we have come to expect 
from time-optimal solutions. As the spacecraft approaches the desired attitude 
the command generated from the error signal approaches zero. Additionally, the 
command generation algorithm does not utilize spacecraft moment of inertia 
information in developing the command. Modern control system engineers have 
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made use of feedforward signals and other more imaginative techniques to 
improve slew performance however; these require an a priori knowledge of the 
maneuver and control solution. This is in effect a step back from the autonomy 
that we seek to achieve. 

2. Singularity Issues 

We have seen that the command generation algorithm outputs a desired 
torque output from the actuator based on the error signals without regard to the 
method of torque generation. This desired torque is sent to the actuator as a 
commanded torque signal to be accomplished. The question to consider is 
under what circumstance is the commanded torque actually output from the 
actuator and what are the effects on the desired maneuver when the 
commanded torque is not output from the actuator. 

Consider an ideal, unlimited torque thruster which is mathematically 
modeled as: 

^"signal, 

”^signal 2 — ~^out 2 
”^signal 3 _ J~out$ 

In this case the input signal (T signal ), results in the corresponding output torque 
{T out ). Therefore, given the desired actuator output torque {T out ), from command 
generation, we can solve for the required input signal to the actuator ( T signal ) as, 

1 0 0 
0 1 0 
0 0 1 

The actuator dynamics, represented in matrix form, in this case as the identity 
matrix, must be inverted to resolve the desired output torque. In this case the 
commanded torque generates the output toque desired. Unfortunately, actuator 
dynamics are seldom this simple. As a result we define an actuator singularity as 
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a condition in which no control torque is generated for the commanded torque 
along a certain direction. Mathematically, this is a condition in which the actuator 
dynamics matrix is singular. 

Consider the case of the asymmetric spacecraft subject to CMG control. 
This case was examined in detail in Chapter VI. Since we have determined that 
the error command generation algorithm is unsuited for generating bang-bang 
controls suppose that we elect to feed forward the bang-bang control structure 
shown in Figure 90. This control torque on the spacecraft would result in a 135 
degree slew about the x-axis for NPSAT 1. 


Inupt Command Torque 



Figure 90 Asymmetric Eigenaxis Rotation Torque Control Solution 

In this case the torque command signal becomes the input to the CMG actuator 
in our spacecraft model (Figure 91). The input or commanded torque signal must 
be resolved into a gimbal angle rate. This is generally accomplished in 
accordance with a command pseudoinverse steering logic 12 . In this steering 
logic for the commanded torque input {u), the CMG momentum rate command 
( h ) is defined as: 
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h = -U —COX/7 

and the gimbal rate command (8 ) is then obtained by:* 

5 =A + h = A T (AA T yh 


CMG Open Loop Command Propagation 



To Workspace 


Figure 91 Spacecraft Model with CMG Actuators 


where the matrix A is a function of the gimbal angles (5 ) and was previously 
defined in Chapter VI (Equation (6.5)). If the matrix A becomes rank deficient for 
certain sets of gimbal angles the pseudoinverse fails to exist and singular states 
are encountered. Previously, we showed that the time-optimal CMG solution 
resulted in well behaved condition numbers for the matrix A (See Figure 88). In 
this case the bang-bang control solution, when used as the desired torque on the 
spacecraft, produced the A matrix condition numbers shown in Figure 92. 
These results are four-orders of magnitude larger than those obtained in the 


* A + is also called the Moore-Penrose inverse of A or the generalized inverse of A. 
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Figure 92 Condition Number of CMG Jacobian Matrix with Independent 

Torque Solution 



Figure 93 Gimbal Angle Rates From Pseudoinverse Steering Logic 

previous section. As a result of the poor condition of this matrix the algorithm 
attempts to drive the gimbal angle rates to infinity. Gimbal angle rates resulting 
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from this simulation are shown in Figure 93. The zoomed in view of the gimbal 
angle rates from 4-6 seconds is shown in Figure 94. 



Figure 94 Gimbal Angle Rates from Pseudoinverse Steering Logic 

(Close-up) 

Gimbal angle rates are driven towards infinity as the singular states are 
approached. Clearly, the singularity is not caused by the bang-bang structure of 
the input torque. Comparing the commanded torque signal (Figure 90) with the 
condition number plot (Figure 92) or the gimbal angle rates (Figure 93) indicates 
that the singular states are encountered prior to the control switch. 

Singularity avoidance is an unintended and yet significant benefit of the 
time-optimal control algorithm. In effect the singularity problem is avoided both 
mathematically and physically by requiring that the time-optimal control vector 
satisfy the actuator constraints in the problem formulation such that, 

u* = -[>4(8 )]S -oox/7(5) 

The actuator dynamics are solved iteratively from left to right vice right to left. 
The previous matrix inversion is eliminated and the solution is optimal with regard 
to the actuator capabilities. 
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In this example we have examined the CMG singularity problem since it is 
well known. However, the same logic could be applied to any actuator system 
that is not an ideal torque generating device. Magnetic torque rods represent a 
second well known system for which the actuator has known singularities. In this 
case, a desired torque output from the torque rod must be resolved in the 
magnetic dipole moment of the torque rod. The torque output of a magnetic 
torque rod is given by: 


f = mxB 
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Therefore, given the desired torque output we have: 
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(7.2) 


The inverse of equation (7.1) does not exist. Therefore it is impossible to 
achieve control of the three body axes by means of the created magnetic dipole 
vector in the satellite 13 . However, using the time-optimal control formulation as 
demonstrated in Chapter III and Chapter IV we require that the optimal control 
vector (m) satisfy equation (7.1). Again, the singular condition is avoided by 
eliminating the matrix inversion and iteratively solving the problem from right-to- 
left. 


C. BELLMAN S PRINCIPLE OF OPTIMALITY 

Bellman’s principle of optimality is a powerful tool when applied to closed 
loop time-optimal control. It states that given an optimal trajectory from a point A 
to a point B, the trajectory to point B from a point C lying on the optimal trajectory 
is also optimal 14 . 
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B 


Figure 95 Principle of Optimality 


In this work this principle serves three purposes. First, it verifies optimality 
in the absence of costate information. Optimal trajectories are obtained and 
verified by recalculating the trajectory from an intermediate point on the 
trajectory. Second, in the presence of a disturbance torque, this principle allows 
the use of previous optimal solutions as guesses for subsequent trajectories. 
Since the disturbed state is near the optimal state using the previously calculated 
state trajectory as a guess for the re-optimized state trajectory significantly 
reduces the computation time. Finally, it allows us to evaluate the accuracy of a 
solution based on a finite number of LGL points. 

Shown in Figure 96 is the now familiar quaternion trajectory for the 
asymmetric spacecraft time-optimal slew maneuver. The intermediate point 
selected is on the optimal trajectory and marked with a bold square. The 
subsequent solution is overlaid with “+” marks. 
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Figure 96 Asymmetric Time-optimal Quaternion Solution - Bellman’s 

Principle of Optimality 

The subsequent solution clearly obeys Bellman’s Principle of Optimality. While 
this may not be sufficient to mathematically demonstrate the optimality of the 
solution it can certainly be considered a necessary condition. The original 
solution was obtained using 30 LGL points (nodes) and was based on a two- 
point guess defined by the initial and final conditions. The subsequent solution 
also used 30 nodes but was provided with a guess which consisted of 15 points 
from the previous solution. The computation timest as determined by the 
MATLAB tic and toe commands are shown in Table 9. 


t All computations were performed on a Pentium 4 Windows based operating system at 3.06 
GHz with 512 MB RAM using MATLAB 6.5.0.1 Release 13. Emphasis is placed on the relative 
vice absolute computational time. 


136 















Computation Time 

Original Solution 

Subsequent Solution 

Reduction, % 

42.5310 seconds 

2.6250 seconds 

93.83 

Table 9 Bellman Principle Efl 

f ect on Computational Time 


The accuracy of the solution obtained is related to the number of nodes 
used in the algorithm. The number of nodes also has a significant effect on 
computational time. Therefore, we strive to choose a number of nodes that 
meets our accuracy goal without creating excessively long computational times. 
As the number of nodes increases, the values returned, states, costates, etc. 
converge so that the difference between solutions approaches zero. Therefore, if 
the number of nodes is increased and the behavior of the solution changes 
considerably the original solution is shown to be inaccurate. In Figure 97, a low 
node (6 nodes) solution is overlaid with a high node solution (30 nodes) from an 
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intermediate point on the optimal trajectory. The significant change visible in the 
state trajectory and final time are indications that the low node solution has 
questionable accuracy. 

D. REAL-TIME OPTIMAL CONTROL 

In this section we provide two examples of closed loop time-optimal 
control. In each case the simulation is as shown in Figure 98. The primary 
tuning parameter for spacecraft performance is the open-loop propagation 
interval. Ignoring computational delay, the state is sampled and an open loop 
control solution is generated time-indexed to the predicted spacecraft state. This 
solution is propagated for a set interval, the propagation interval, after which the 
state is re-sampled and a re-optimized trajectory is computed. Currently, 
computation times are recorded for comparison and evaluation but ignored in the 
simulation. 



Figure 98 Closed Loop Simulation Model 


1. Asymmetric Spacecraft with Idealized Actuator Control 

The time-optimal reorientation of an asymmetric spacecraft with idealized 
actuator control was addressed in detail in Chapter III. The spacecraft and 
orbital parameters remain unchanged for the closed loop simulation and are 
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available in Chapter III. A constant disturbance torque has been added to the 
model. This disturbance torque will perturb the model off the optimal trajectory 
and require re-optimization. Since we know from the previous chapter that the 
time-optimal open loop maneuver is approximately 6 seconds the disturbance 
torque has been amplified from what might be considered normal for graphical 
visualization. The disturbance torque has been selected as: 

M d =0.108 N.m 

M d =0.108 N*m 

u x 

M d =0.0317 N*m 

u x 

The open loop time-optimal control solution and the propagated solution in 
the presence of this disturbance torque are shown in Figure 99. The solid lines 
indicate the optimal trajectory; the propagated trajectory is marked as “+’. It is 


Quaternion History Solution & Propagation Angular Rate History Solution & Propagation 



Figure 99 Asymmetric Spacecraft Time-optimal State Solution and 

Propagation 

clear that in presence of a disturbance torque the open loop solution would not 
deliver acceptable performance. 

The propagation interval for this closed loop simulation was selected as 1 
second. The results of the first interval of propagation are shown in Figure 100 
and Figure 101. At this time the spacecraft state is sampled and the time-optimal 
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Figure 100 Asymmetric Spacecraft Quaternion Propagation Interval 1 



Figure 101 Asymmetric Spacecraft Angular Rate Propagation Interval 1 


solution is computed from the current attitude and rate to the desired attitude and 

rate. In the next two figures, figure 102 and Figure 103) the re-optimized 

trajectory is shown along with the original trajectory. The basic structure of the 

optimal state trajectories is the same and the end point conditions are met. 
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Figure 102 Asymmetric Spacecraft Quaternion Propagation Interval 2 



Figure 103 Asymmetric Spacecraft Angular Rate Propagation Interval 2 

The sample and update cycle continues from Figure 104 through Figure 
113. Each time the state is sampled the time-optimal solution is generated and 
propagation continues with the updated solution. 
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Figure 104 Asymmetric Spacecraft Quaternion Propagation Interval 3 



Figure 105 Asymmetric Space craft Angular Rate Propagation Interval 3 


With each update it is clear that a re-optimized trajectory is generated from the 
current state of the spacecraft to the final state which remains a constant. 
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Figure 106 Asymmetric Spacecraft Quaternion Propagation Interval 4 



Figure 107 Asymmetric Spacecraft Angular Rate Propagation Interval 4 
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Quaternion Solution and Propagation 



Figure 108 Asymmetric Spacecraft Quaternion Propagation Interval 5 



Figure 109 Asymmetric Spacecraft Angular Rate Propagation Interval 5 
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Figure 110 Asymmetric Spacecraft Quaternion Propagation Interval 6 



Figure 111 Asymmetric Spacecraft Angular Rate Propagation Interval 6 
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Figure 112 Asymmetric Spacecraft Quaternion Propagation Interval 7 



Figure 113 Asymmetric Spacecraft Angular Rate Propagation Interval 7 


The completed maneuver is shown in Figure 112 and Figure 113. The 
time required to complete the 135 degree reorientation about the x-axis in the 
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presence of the disturbance torque is 7.2652 seconds. The original solution 
time, modeled without a disturbance was 6.1825 seconds. Each subsequent 
solution calculated from the sampled state met the feasibility and optimality 
criteria that were established for this problem formulation in Chapter III. 

Computation time is related to propagation interval in that longer 
propagation intervals in the disturbance field cause the spacecraft to deviate 
further from the optimal trajectory. This results in longer computation times. The 
computation times for the simulation are shown in Figure 114. The average 
value of the update solutions (solutions 2 through 7) is 11.9296 seconds. This is 
due to the large disturbance torque modeled to aid in the visual clarity of the 
graphics. When the disturbance torque is modeled as something more realistic 
for the spacecraft and orbital parameters the previous disturbances are reduced 
by 5 orders of magnitude. The resulting update computation times then 
decrease to an average value of 4.0566 seconds. 



Figure 114 Asymmetric Spacecraft Idealized Actuator Time-optimal 

Solution Computation Time 
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a. Numerical Considerations and Notes 

Simulated spacecraft states were numerically generated by a 
MATLAB ODE45 propagation subroutine. The generated control vector plus the 
disturbance torque were used as input; the spacecraft state was extracted 
following propagation to simulate a sensor reading. In some cases due to 
numerical approximations and truncation errors the norm of the quaternion was 
not equal to one to the precision required for the optimization algorithm (DIDO) to 
determine a feasible solution. In order to resolve this problem the sampled-state 
quaternion vector was normalized as follows: 


' nomai |<4 

Additionally, Figure 114 indicates, and it was observed, that the 
computation time for the final re-optimization of a trajectory was in general longer 
than previous computation times. Since this algorithm was originally conceived 
for large angle slew maneuvers constant scaling parameters were applied based 
on the maneuver open loop solution. As the maneuver approaches completion 
and the trajectory is re-optimized the magnitude of the problem has changed to 
the extent that the original scaling parameters are poorly chosen. The problem 
requires a scaling parameter adjustment based on the remaining horizon. 

2. Asymmetric Spacecraft with Magnetic Torque Rod Control 

The time-optimal reorientation of an asymmetric spacecraft under 
magnetic torque control was treated in detail in Chapter IV. Once again the 
spacecraft and orbital parameters are based on the NPSAT 1 small satellite 
program and found in Table 5. 

The closed loop simulation, while made more complex by the time-varying 
magnetic field, is accomplished as shown in Figure 98. The constant disturbance 
torque selected for propagation is: 
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M d = 5.4x1 CT 6 N.m 

u x 

M d =5.508x10 6 N*m 

y 

M d =6.34x10 7 N.m 

dz 

The larger disturbance torque, used in the previous example for visual clarity, 
would overpower the torque available from the system and make the spacecraft 
uncontrollable. Recall that the time-optimal open-loop maneuver required 
approximately 273 seconds. The state trajectories resulting from the propagation 
of the open loop solution in the presence of the disturbance torque are shown in 
Figure 115 and Figure 116. The end point error, although difficult to see in the 
graphical presentation is significant enough to warrant correction. A comparison 
of the actual vs. desired state is shown in Table 10. 



Figure 115 Magnetic Torque Open Loop Solution Propagated with 

Disturbance Torque 
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Figure 116 Magnetic Torque Open Loop Solution Propagated with 

Disturbance Torque (2) 


Propagation Error with Disturbance 

State 

Desired 

Actual 

Error 

ql 

0.9239 

0.9216 

0.0023 

q2 

0 

-0.01 

0.01 

q3 

0 

0.0206 

0.0206 

q4 

0.3827 

0.3879 

0.0052 

Omega 1 

0.00 

-1.00E-04 

1.00E-04 

Omega 2 

7.72E-04 

7.00E-04 

7.25E-05 

Omega 3 

7.72E-04 

1.30E-03 

5.28E-04 


Table 10 State Trajectory Propagation Error with Disturbance 


A propagation interval of 15 seconds was selected for this simulation 
based on the open loop maneuver time and the magnitude of the disturbance 
torque. As before, following propagation the spacecraft state is sampled and the 
trajectory is re-optimized from the current state. 
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DIDO has internal constraints which define the accuracy with which the 
optimal trajectory must reach the defined end points. In order to facilitate the 
solution, box constraints were applied to the desired end state. These form an 
upper and lower bound on the end state such that and acceptable end state is ± 
epsilon from the exact end state. For this simulation epsilon was set as, 

q E =0.006 
co E =0.0001 

The end state box constraints appeared to increase computational speed. 
In Figure 117, the first update to the disturbed quaternion trajectory is shown and 
appears to complete the maneuver faster than the original undisturbed solution. 
This is not surprising. The open loop solution was formulated to the end point 
and the subsequent solution was formulated to the end point with the addition of 
the box constraint. Additionally, it is possible that the constant disturbance 
torque is initially providing some benefit to maneuver speed. Flowever, the final 
closed loop state trajectories figure 118 and Figure 119) demonstrate that to 
meet the established accuracy requirements in the presence of the disturbance 
torque requires additional time. 
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Quaternion Solution and Propagation 



Figure 117 Magnetic Torque Quaternion Propagation Interval 2 



Figure 118 Magnetic Torque Quaternion Propagation - Maneuver 

Completion 
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Figure 119 Magnetic Torque Angular Rate Propagation - Maneuver 

Completion 

A comparison of the open loop maneuver with the disturbed closed loop 
maneuver shows that the disturbance increases the overall maneuver time. 

The disturbance free simulation required 273.4117 seconds to complete the 135 
degree x-axis slew. The close loop maneuver in the presence of the disturbance 
torque required 285.3174 seconds. The overall accuracy improvement is shown 
in Table 11. The improvement appears significant enough to warrant the 
application of closed loop control. Again, all subsequent solution updates met 
the feasibility and optimality necessary conditions that were established in 
Chapter IV. 


Closed Loop Propagation Error 

State 

Desired 

Actual 

Error 

Error Reduction, % 

ql 

0.9239 

0.9263 

0.0024 

4.35 

q2 

0 

0.0042 

0.0042 

58.00 

q3 

0 

0.006 

0.006 

70.87 

q4 

0.3827 

0.3767 

0.006 

15.38 

Omega 1 

0.00 

-1.00E-04 

1.00E-04 

0.00 

Omega 2 

7.72E-04 

8.00E-04 

2.75E-05 

62.01 

Omega 3 

7.72E-04 

7.00E-04 

7.25E-05 

86.26 


Table 11 Closed Loop State Propagation Results 
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As before, computation time is the key to establishing the feasibility of 
closed loop real-time optimal control. The computation of magnetic torque time- 
optimal maneuvers is significantly complicated by the time-varying magnetic 
fields. This is illustrated by the long computation time required for the open loop 
solution shown in Figure 120. The open loop solution was based on a two-point 
guess, the end points, and calculated over 30 nodes. The subsequent solutions 



Figure 120 Magnetic Torque Time-optimal Maneuver Computation Time 

were also calculated over 30 nodes but used the previous state, control solution 
as a guess. The time required for the open loop solution was 324.9680 seconds, 
subsequent solutions had an average time of 12.3407 seconds. This represents 
a computation time reduction of 96.2 percent. 


E. CONCLUSIONS 

This section has demonstrated that classical feedback control theory is ill 
equipped to handle time-optimal maneuvers. More recent techniques employing 
neighboring optimal control laws have improved slew performance but still 
require inverting the actuator dynamics and therefore require singularity 
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avoidance. By requiring that the control solution satisfy the actuator dynamic 
constraint equations singularity avoidance is automatically eliminated. 

Strizzi et al. 15 noted, and engineering common sense agrees, that if open 
loop control solutions can be generated in real-time they become equivalent to 
feedback control laws. Computational power and numerical methods have come 
together to make this a reality. 

The two numerical examples presented in this section illustrate the 
concept of fast and slow dynamics. In the idealized actuator example, a 
computational time of 10 seconds would not suffice for real-time control. 
However, that computation time applied to the magnetic torque example would 
seem an excessive requirement. The computational speed required to achieve 
real-time optimal control is a product of both the plant dynamics and desired 
performance characteristics. 

Empirically it was observed that the time required to calculate the 
subsequent trajectories is a function of two factors, deviation from optimal 
trajectory and amount maneuver remaining. As the frequency at which updates 
are made increases, the amount of deviation from the optimal trajectory due to 
disturbance torques and plant uncertainties necessarily decreases. Therefore, 
computation time required decreases as update frequency increases. Second, 
as the maneuver approaches completion the computation time required 
asymptotically approaches a minimum value. 

Using previous solutions to jump start the updated solution was shown to 
significantly reduce solution computation time. Further reductions in actual 
computation time are predicted by reducing operating system overhead and 
employing field programmable gate arrays. 
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VIII. FUTURE WORK 


The work in this thesis is by no means ready for on-orbit testing. What 
follows are areas of improvement and future research that time precluded. 


A. ALGORITHM IMPROVEMENT AND OPTIMIZATION 

This algorithm was originally conceived for large angle slew maneuvers. 
To that end, numerical scaling was used to adjust all parameters to within the 
same order of magnitude. This improved the accuracy of results and reduced 
computational time. However, as the spacecraft moment of inertia was varied 
there were conditions under which separate scaling for the angular rates of the 
body axes would have improved performance. A scaling format based on 
structure array variables 1 would simplify coding and allow separate scaling for 
each variable. 

Additionally, the time scaling was designed for the initial open loop 
solution and held constant throughout the closed loop simulation. Therefore, as 
the time-to-go horizon decreased, the problem was poorly scaled in some cases. 
As the spacecraft state is sampled for each re-optimization the time scale should 
be adjusted so as to ensure that subsequent solutions are well scaled for 
accuracy and to reduce computational time. 

Though great emphasis has been placed on computational time in the 
realization of real-time optimal control, the algorithms used in this work were not 
designed to be time-efficient. Some reconfiguration of the loops and control 
structure of the algorithms would undoubtedly result in additional computational 
speed. 

Finally, though not the subject of this work, increasing the speed of the 
reusable software package DIDO is being investigated. Through encoding 
analytic Jacobian significant computational speed enhancements are possible . 2 - 3 
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B. COMPUTATIONAL DELAY MODELING 


Key to the implementation of this algorithm is the computational delay 
associated with the solution to the optimal control problem. As shown in Figure 
121 the computational delay (S c ) is the time from when the state is observed at ti 

until the time when the re-optimized control solution is available at h. During that 
time the state is continuing to propagate under the previous control solution. 



Figure 121 Conceptual Computational Delay Model 
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Bellman’s principle of optimality says that in the absence of a disturbance 
torque the sampled state at ti would lie on the optimal trajectory and the updated 
control solution would be identical to the original solution. However, in the 
presence of a disturbance, the sampled state could be off the optimal trajectory 
and the updated control solution would not be applied until some non-zero time 
later. In this case, the re-optimized solution would be applied to a system whose 
state has propagated to another point, presumably off the re-optimized trajectory 
thus inducing some error into the system. In this work, the computational delay 
was taken to be negligibly small. In fact, the effect of the computational delay is 
a function of both the spacecraft dynamics and the desired performance. 
Numerical modeling and simulation of the effect of computational delay is 
necessary to determine both a sample rate and maximum acceptable 
computational delay. 

C. EXTENSION TO HARDWARE 

Finally, much research is needed, and has begun, to implement these 
algorithms in hardware. Implementation via field programmable gate arrays 
(FPGA) is in ongoing. Through the efforts of others this project will see air¬ 
bearing testing and eventual on-orbit testing. 
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